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NUMERICAL APPROXIMATION OF BOUNDARY CONDITIONS WITH 


APPLICATIONS TO INVISCID EQUATIONS OF GAS DYNAMICS 

H. C. Yee 

Ames Research Center 


SUMMARY 


A comprehensive overview of the state of the art of well-posedness and 
stability analysis of difference approximations for initial boundary value 
problems of the hyperbolic type is presented. The applicability of recent 
theoretical developments to practical calculations for nonlinear gas dynamics 
is examined. The one-dimensional inviscid gas-dynamics equations in 
conservation-lav-form are selected for numerical experiments. The class of 
implicit schemes developed from linear multistep methods in ordinary differ- 
ential equations is chosen and the use of linear extrapolation as an explicit 
boundary scheme is emphasized. Specification of boundary data in the primitive 
variables and computation in terms of the conservative variables in the 
interior are discussed. Some numerical examples for the quasi-one-dimensional 
nozzle are given. 


1 . INTRODUCTION 


The proper specification of boundary conditions which yield a well-posed 
problem for a partial differential equation is essential for the behavior of 
the solution. Overspecification of boundary data precludes the existence of 
smooth solutions except in very special unrealistic situations in which the 
exact solution is known on the boundary without error. In the development of 
difference approximations for mixed initial boundary value problems in the 
applied science field, the boundary conditions may be quite difficult to con- 
struct, and a poor choice can lead to inaccuracies and instabilities. Part of 
the difficulty starts with the original differential equations where the 
proper boundary conditions are not always known (nonlinear fluid dynamics 
problems, for example). The problem is compounded in the difference schemes 
where quite often extra boundary conditions are needed because the difference 
equations are of higher order than the differential equations. Therefore, in 
the study of how the extra boundary conditions affect the stability and accu- 
racy of numerical schemes, we not only have to examine the difference schemes 
used, but we also have to first examine the well-posedness of the original 
differential equations (refs. 1-15). Thus a good understanding of the theory 
of "well-posed problems" is a necessity. 

The two principal objectives of this report are to (1) present a compre- 
hensive overview of the state of the art of well-posedness and of stability 
analysis c: difference approximations for initial boundary value problems of 
the hvperbolic type, and (2) to examine the applicability of the current theorv 



to the inviscid (Euler) equations of gas dynamics. (We will us he terms 
”inviscid gas-dynamics equations” and ”Euler equations of gas d *tmics" inter- 
changeably.) Through an understanding of the theory, we can g some insights 
into how to impose the physical boundary conditions more corrc /, and we can 
be guided in the construction of stable numerical schemes for practical 
problems. In this context, ”stable numerical schemes” are sc les that are 
stable for the combined interior and boundary schemes. Read s who are famil- 
iar with the theory and who are only interested in the appli cion can skip 
the first four subsections of the second section and can sk:^ the third section 
altogether. 

In this report, we will discuss several ways of formulating the boundary 
approximation for the one-dimensional inviscid gas-dynamics equations in con- 
servative form. Since in general the Euler equations have mixed positive and 
negative eigenvalues , appropriate one-sided and uncentered boundary approxima- 
tions are essential. Some of the methods proposed in this report combine the 
theory of Gustafsson et al. (ref. 9) with the flux-vector splitting technique 
of Steger and Wanning (ref. 16) to study the applicability of some uncondi- 
tionally stable schemes for the one-dimensional (1-D) linearized Euler equa- 
tions to their nonlinear counterpart- A few detailed numerical results for 
the quasi-l-D nozzle with various inflow-outflow conditions are given. The 
boundary approximations being used are one-sided spatial differencing and 
linear extrapolation. It was found that we can use fairly large CFL numbers 
(i.e., Courant, Fredrick, and Levy condition for the stability of differences 
schemes) . 

The review of the theory of well-posed problems and stability analysis of 
difference schemes is desirable because significant progress on a general, 
workable theory for the initial boundary value problem of the hyperbolic (and 
parabolic) type is quite recent. Much of it begins with the work of Kreiss 
(ref. 1) published in 1970. The recent research papers on this rapidly- 
developing subject are principally addressed to highly-theoretical audiences, 
and there is no text or basic, up-to-date review article covering the material. 
A primary purpose of this report is to collect the relevant information and to 
identify some of the strengths and weaknesses of the existing theory when it 
is applied to physical problems. The material is presented with the needs of 
applied scientists in computational fields in mind. Consequently, basic con- 
cepts and practical ideas are emphasized while exact mathematical definitions 
and theorems are minimized. Only initial boundary value problems of the 
hyperbolic type are considered. 

Section two of this report is a review of the state of the art of how to 
impose boundary conditions in order to obtain a well-posed problem. Section 
three is a cdmpfehehsive review of the current status of stability analysis of 
difference approximations. For example, a recent result by Oliger (ref. 10) 
provides a useful guide in the construction of composite stable schemes. In 
the fourth section, a detailed application of these theories is given for the 
one-dimensional Euler equations of gas dynamics; several numerical experiments 
are included* In addition to the numbered references that are cited through- 
out the text, a separate bibliography is provided. The bibliographic entries 
are categorized according to their particular relevance to sections 2, 3, 
and 4. 
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2. WELL-POSEDNESS OF INITIAL BOUNDARY VALUE PROBLEMS (IBVP) 
FOR HYPERBOLIC EQUATIONS 


In this section, we summarize the status of well-posedness of initial 
boundary value problems (IBVP) for hyperbolic partial differential equations. 

A more detailed discussion of this subject is given in subsequent sections. 

Some of the related mathematical definitions and examples are described in 
appendixes A and B. Readers who are not familiar with the definitions should 
refer to appendix A. 

The term "well-posed," or correctly-posed, problem appears frequently in 
the literature. There are many different definitions for a well-posed partial 
differential equation; for example, well-posedness in Hadaraand’s sense is dif- 
ferent from well-posedness in Kreiss’s or Petrovskii’s sense. The various 
definitions can be found in sources contained in the first two sections of the 
Bibliography. In this report, we only consider well-posedness in Kreiss’s 
sense; that is, "well-posedness in the L 2 norm." The basic requirement for 
a well-posed IBVP is to not overspecify or underspecify the boundary conditions 
with given smooth initial data. In order to mathematically define a well- 
posed IBVP, we have to establish the existence and uniqueness of the solution 
and its continuous dependence on the initial and boundary data or to establish 
the existence of certain a priori estimates or energy inequalities, 

Well-posedness of the governing partial differential equation is a very 
crucial consideration commonly overlooked by investigators in the field of 
computations; that is, the problem is defined only when a proper set of initial 
and/or boundary conditions is given. We cannot expect our difference approxi- 
mations to be reasonable if they approximate a problem that does not have 
reasonable solutions. In many instances, a good understanding of the theory 
of well-posed problems can guide us to exclude many boundary conditions which 
might look physically reasonable. 

The theory for the IBVP of 1-D systems or degenerate 1-D systems (higher 
dimension systems that can be reduced to 1-D problems (ref. 17)) has been 
established for some time. For higher dimension systems (with constant coef- 
ficient problems) , results are known for the strictly hyperbolic and the sym- 
metric hyperbolic case (see "More Than One Space Dimension," p, 11, for defini- 
tion) . Some partial results for the multidimensional Euler equations were 
established by Oliger and Sundstrdm (ref. 6). 

The following sections are summary discussions of ways to impose or to 
check for well-posedness of IBVP for hyperbolic equations in the L 2 norm 
(see appendix A). We will discuss the following types of problems; 
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1 . 1-D scalar equation 

2. 1-D system of equations 

3. Several space dimensions equations 

with constant coefficients, variable coefficients, and quasilinear proper s. 

We assume that the problems we are considering have smooch initial data, 1 

of the permissible ways of imposing boundary conditions in the subsequer 
sections are necessary and sufficient conditions for well-posedness of 1-D 

hyperbolic equations. The discussions are based on the method of chare ris- 

tics. For the more-than-one-space-dimensions problem, the analogous fc ula- 
tion need not be well-posed. A proper way of getting a necessary and suffi- 
cient condition in this case is by the normal mode analysis (ref. 1). One way 
of getting a sufficient condition is by the energy method (refs. 4, 12). 


Scalar Equation 

Consider the problem 

Uj, + cu^ = 0 t > 0 , c real constant (la) 

with initial condition 

u(x,0) - f(x) (lb) 

We can divide the above problem into the following three categories. 

a. The Cauchy (initial value) problem (-°° < x < <») ; The exact solution 
is given by 

u(x,t) = f(x - ct) (2) 

Hence the solution of (1) is constant along the characteristic lines 
X - ct = constant. There are no boundary conditions involved since -oo<x<®. 

b. Half-space problem (0 < x < °°) : 
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u (x, 0) = f(x) 
c > 0 

Figure 1.- Half-space problem (c > 0) . 



Figure 2.- Half-space problem (c < 0) . 


If c > 0, then u(x,t) is only determined in the triangular region 
X - ct > 0 (see fig. 1). In this case, we need a boundary condition 

u(x,0) » g(t) t > 0 

to determine the solution for x - ct < 0. If c < 0, then u(x,t) is 
uniquely determined by (2) and it is not appropriate to specify a boundary 
condition at x “ 0 (see fig. 2). Note that the solution u(x,t) is continu- 
ous in a neighborhood of x - ct “ 0 if and only if f and g are continuous 
and satisfy the compatibility condition 

f(0) - g(0) 

c. Finite domain problem (0 < x < 1): 



0 1 


u(x, 0) = f(x) 
c > 0 

Figure 3.- Finite domain problem 
(c > 0). 



u(x, 0) = f(x) 
c< 0 

Figure 4.- Finite domain problem 
(c < 0). 
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In this case (see figs, 3 and 4), the necessary and sufficient boundary con- 
dition to produce a well-posed problem is 


u(0,t) » g(t) 

if 

c > 0 

u(l,t) - g(t) 

if 

c < 0 


System of Equations 

A system of first-order constant coefficient partial differential 
equations 

ut + Au^ "0 t>0, 0<x<l 

is said to be hyperbolic if A is diagonalizable and with real eigenvalues. 
We can divide the system of equations into the following five categories; 

a. System of hyperbolic equations in diagonal form 

b. System of hyperbolic equations in coupled form 

c. Nonpositive definite systems 

d. Variable coefficients 

e. Quasi-linear systems 

Each is discussed below. 


System of hyperbolic equations in diagonal form : 

^t t:>0, 0<x<l 

or 



(3) 
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are positive definite diagonal matrices. We can categorize the system further. 

i) N decoupled equations with decoupled boundary conditions: The solu- 

tion is uniquely determined if we specify initial values 

u^(x,0) « f^(x) 

u^^(x.O) - f^^(x) 

and boundary conditions 

u^(0,t) » g^(t) 
u^^(l,t) - g^^(t) 

In this case, we just solve N independent scalar equations. 

ii) N decoupled equations with coupled boundary conditions: We can 

couple these equations if we replace the above boundary condition by 

u^(O.t) = SjU^^(O.t) + g^(t) 

t > 0 . (4) 

u (l,t) = Sj.j.u^(l,t) + g^^(t) 

where Sj. Sjj are rectangular matrices with dimension £ x (n - 5.) and 
(N - £) X ji, respectively. From the examination of how the direction of the 
characteristic lines (and the initial data) determine the solution of the 
finite domain scalar equation, we can conclude that the solution of equa- 
tion (3) is again uniquely determined by conditions (4) and the initial data. 
Geometrically the values of u^ are transported along the characteristic to 
the boundary x = 1 (see 5). Then, by the boundary conditions 

= Snu-‘-(l,t) + g-‘-^(t), these values are transformed into values for 
u , which are then transported to the boundary x = 0, etc. Therefore, the 
number of boundary conditions for x « 0 is equal to the number of positive 
eigenvalues of the matrix A. And the number of boundary conditions for 
X - 1 is equal to the number of negative eigenvalues. Thus a necessary and 
sufficient condition for the IBVP of system (3) to be well-posed is to impose 
the boundary conditions in the form of (4). But the analogous formulation for 
problems in more than one space dimension is not necessarily well-posed. This 
subject is discussed briefly in the next subsection. 
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0 1 
u(x,0| =f(x) 


Figure 5.-- Coupled boundary conditions. 


System of hyperbolic equations of coupled form : In most applications, the 
system of differential equations is coupled 

^t ^ t>0, 0<x<l 

where A is assumed to be a constant matrix which can be diagonalized by a 
transformation matrix T. 


T"^AT = A 

where A has the same form as (3) (or we can rearrange A In order to have 
the same form as (3)). For x ■ 0, the boundary conditions consist of I 
lln-jar relations among the components of u, that Is, In matrix form 

Lu(0,t) =■ g(t) (5) 

where L Is an J, x N matrix. Recall that I Is the number of positive 
eigenvalues of A. Let the characteristic variables be defined by 

w - T“^u 

Then w is the solution of equation (3) , and the problem Is well-posed if the 
boundary condition 


LTw(0,t) - g(t) 

at X = 0, and the similar boundary condition at x - 1, can be written in the 
form (4). Here, the rank of L must be equal to the number of positive 
eigenvalues of A at x * 0. Therefore, any boundary conditions specified 
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for the original system, must be transformable to boundary conditions of the 
form (4). 

Nonpositive definite systems ; If or are not positive definite, 

then the components uJ(x,t) corresponding to Xj - 0 must be considered as 
outgoing variables and will be included In for x = 0 and in u^ for 

X » 1. (Variables associated with negative eigenvalues are termed as outgoing 
variables, and variables associated with positive eigenvalues are termed as 
incoming variables for the left boundary.) Since the characteristic is verti- 
cal, the solution along this characteristic is determined by the initial con- 
dition, Therefore, our discussion can assume that A^ > 0 for x - 0 and 
A >0 for x “ 1. That is, we should not specify the corresponding jth 
boundary condition with respect to Xj ■ 0 . 

An example of a well-posed system of hyperbolic equations (from Kreiss 
and Oliger, ref. 2) follows. Consider 



The eigenvalues, Xj of A, are 

Xj - -c 

Xj - -(C + 1) 

X 3 - -(c - 1 ) 

Assume 0 < c < 1, then A has one positive and two negative eigenvalues. 
Therefore, we have to specify one boundary condition at x ■ 0 and two 
boundary conditions at x » 1 . 

Let us consider the boundary conditions 


U3(0.t) - 0 ' 

Ui(l.t) - 0 . (6) 

ujd.t) - 0 j 

and check whether ( 6 ) produces a well-posed problem. That is, we have to show 
that these conditions are, after transformation, of the form (4). The trans- 
formation T that diagonalizes A is 
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and (6) becomes 


WjCO.t) - W2(0,t) 
w^(l,t) - 0 
w^Cl.t) » WgCl^t) 


Here 


and 


■ W 3 




0 


The above equations are obviously of the form of equation (4). 


Variable coefficients : 


+ A(x,t)Ujj - 0 t>0, 0<x<l 

For every fixed t « tg , the form of the well-posed boundary conditions at the 
boundaries x * 0 and x ■ 1 is determined by the systems with constant 
coefficients 

Vj. + A(0,tg)Vj^ - 0 

Wj. + A(l,tQ)w^ - 0 

respectively. This is the so-called "freezing" method. Note that the theory 
does not cover the case when an eigenvalue of A(0,t) or A(l,t) changes sign 
in the time interval of interest. Therefore, any eigenvalue of A(0,t) or 
A(l,t) has to remain the same sign over the time interval of interest if the 
current theory is to be applied, 

Quasilinear systems ; 

Uj. + A(u,x,t)u^ - 0 (7) 

Assume A(u,x,t) is a smooth function of the arguments u, x, and t, and that 
u can be represented as 
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u(x,t) * U(x,t) + u(x,t) 


where U(x,t) represents a known smooth solution and u(x,t) a small perturba- 
tion. Linearizing equation (7) with respect to u(x,t) gives us a linear 
system 

u^ + A(U,x,t)u^ - B(U,x,t)u + F(U,x,t) 

The boundary conditions for the well-posed problem are determined by A(U,x,t) 
which is discussed in the variable coefficient case. The matrix B(U,x,t) only 
affects the initial conditions and F is the nonhomogeneous part of the 
equations . 


More Than One Space Dimension 
Half-space scalar problems (x > Q» -” < y < «) ; 

+ au^ + bUy *0 x>0, -“<y<®,t> 

u(x,y,0) = f(x,y) a, b real 
The solution of (8) is 

u(x,y,t) - f(x - at, y - bt) 

If a < 0, then u(x,y,t) is completely determined by f. If 
have to specify boundary values 

u(0,y,t) » g(y»t) 

Again, the sign of "a" determines whether we need to impose boundary values. 

The additional space dimension y does not interfere with the above boundary 
condition, since the y domain is -« < y < «. 

Bounded region (scalar problem) : Consider equation (8) in a closed bounded 
region U with smooth boundary 


( 8 ) 


a > 0 then we 
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Let (x*,y*) be a point on The boundary d^ta which should be specified at 

(x*,y*) are again determined by the direction of the characteristic at t 
point. That is, we have to pointwise map the boundary onto the tangen ^ane 

(n,T). This can be formalized by introducing a new coordinate system 
origin (x*,y*) and axes directed as the tangent t and the internal r al n 

X * (x - x*)cos 0 - (y - y*)sin 0 

y - (x - x*)sin 0 + (y - y*)cos 9 

where 0 is the angle between the y and t axis. The new transfer ed equa- 
tion has the form 

+ au- + bu- - 0 5t>0, -«<5^<a>,t>0 


where 


a ■ a cos 0 - b sin 9 
b ■ a sin 0 + b cos 0 

The sign of a determines whether we need to specify boundary data at (x*,y*). 

Mo re-than-one- 3 pace-dimens ion system of equations : The form of the 

necessary and sufficient conditions for well-posed problems for the 1-D system 
only give necessary conditions for their multidimension counterpart. In order 
to obtain necessary and sufficient conditions; we have to resort to normal- 
mode analysis (see appendix A for definition) and Laplace transform (refs. 1-3) 
types of approach. Known theory by the normal mode analysis is only for 
strictly hyperbolic systems and symmetric hyperbolic systems. 

Consider a first-order system in two space dimensions 

u^ + Au^ + Buy ■ 0 X > 0, -« < y < ®, t > 0 

with constant coefficient matrices A and B with dimension N x N. The sys- 
tem is hyperbolic if for all real u >2 with - 1, there is a 

nonsingular transformation T * T(u)i, 0 D 2 ) for which both T and T“^ are 
uniformly bounded such that 


T(u>^A + uj2B)T“^ - 


^2 


0 



and Xj are real. If all the Xj are distinct, the system is strictly 
hyperbolic. If the matrices A and B are both symmetric, then we call the 
system symmetric hyperbolic. We remark that the 2-D and 3-D inviscid gas- 
dynamics equations are not strictly hyperbolic. 
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We will not discuss the normal mode analysis method here; interested 
readers are referred to Kreiss’s original paper. Here, we want to discuss the 
example that Kreiss and Oliger (ref, 2) and Kreiss (ref. 3) have used to 
illustrate the insufficiency of the method of characteristics, Kreiss has 
considered the linearized shallow-water equations 

w^ + Aw^ -I- BWy = 0 x>0, -°°<y<o°,t>0 

where 



Then the matrix A has two positive eigenvalues and two boundary conditions 
have to be described at x = 0. Kreiss used boundary conditions of the form 

V = 0 

Bu + a(^ = 0 

Choosing different values of a,B, the system can have solutions (1) that grow 
arbitrarily fast with time, (2) that have too much reflection at the boundary, 
or (3) that are smooth and well-behaved. The following are his findings: 


i) 

For 

a 

< -1, S = 1 

, situation (1) 

occurs 

ii) 

For 

B 

= 0, a = 1 , 

situation (2) 

occurs 

iii) 

For 

B 

o 

II 

II 

and a = B = 1 

situation (3) occurs 


For problems in several space dimensions that have smooth boundaries and 
smooth coefficients, Majda and Osher (ref. 5) showed that we only need to look 
at the family of frozen constant-coefficient problems on half-spaces obtained 
by locally mapping the boundary onto the tangent plane at each point of the 
boundary, freezing the coefficients locally and disregarding the rest of the 
boundary. They showed that the original problem is well-posed if every member 
of this family of problems is well-posed. 
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Gas -Dynamics Problems 


For the inviscid systems with smooth solutions or problems with low 
Reynolds number, Oliger and SunstrOm (ref* 6) and Oliger (ref* 7) have est 
lished conditions for well-posedness of multidimensional problems. For s 
sonic inflow problems this set of admissible conditions, with a few excel 
is of the form similar to (5) with almost full nonzero entries for L. ^ t 
means we have to impose a set of conditions that are linear combinations the 

physical variables instead of the physical variables themselves. But, i 
physical reasons, we can only specify boundary data that are measurable. In 
this case, most of the admissible boundary conditions for subsonic inflow 
problems do not have physical significance or are not measurable quantities. 

For example, specifying pressure for subsonic inflow is very desirable physi- 
cally, but theoretically the solution results in continued loss of smoothness 
globally. The most physically well-posed boundary conditions they have shown 
are when all of the velocity components along with either the density or the 
temperature are specified. By using the method of characteristics or the 
normal mode analysis, the 1-D inviscid gas-dynamics equations possess some 
features that their higher dimensional counterparts do not have; that is, 
there are boundary conditions that are well-posed for the 1-D Euler equation 
but are not well-posed for the 2-D and 3-D case. This 1-D case is discussed 
in more detail in the next section. 


One-Dimensional Inviscid Equations of Gas Dynamics 


In one spatial dimension, the inviscid equations of gas dynamics can be 
written in the conservative form as 


where 


au ap(u) . . 

at ax ■ ^ 



are the conservative variables and 


F 


m 


m^/p -h p 
i(e -f p)m/pi 


(9) 


is the flux vector, and m “ pu. The primitive variables (denoted by U) are 
the density p, the velocity u, and the pressure p. The total energy per 
unit volume, e, is defined as 
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e = pe + pu^/2 

with e as the internal energy per unit mass* The pressure p for a perfect 
gas is defined as 

p = (y - 1) [e - mV2p] 

where y is the ratio of specific heats. We can write equation (9) in quasi- 
linear nonconservative form as 


0 

at 3x 



where 

A = M’^AM 



9x 3x 
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We can freeze system (11) (assume constant values of A ■ Aq) (notice t. . we 
do not have to freeze the coefficient before getting into this form; th 
freezing of the coefficient is for later analysis) and transform (11) 


T-i 



^A^TT-^U^ 


0 


where 


and 



with W as the characteristic variables 



W - T"^U 


( 12 ) 


and Uq, Cq, and are the "frozen coefficient" values or, numerically, the 
values at a given time-step and grid point. System (12) is transformed into 
the following characteristic form 



On the other hand, we can locally linearize system (11) into 

1^ + A(Uq) 1^ + B(Uq)U + F(Uj,,x,t) = 0 - (13) 

where U * + U, and represents a smooth solution and U is a small 

perturbation. This local linearization of (11) is for checking the well- 
posedness of boundary conditions. The boundary conditions are then determined 
by A[Uo(x,t)]; that is, the form of the boundary conditions at the boundaries, 
say, for fixed t-t^, 0<x<l, are determined by the systems with constant 
coefficients 

1^ + A[Uo(0,tQ)] - 0 at x-0 (14) 
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3w 

at 


at 


X = 1 


(15) 


+ A[D„a.t„)l-U-0 

But if we know the type of inflow-outflow conditions beforehand, there is a 
very simple way of checking the well-posedness (instead of eqs. (14) and (15)) 
once a given set of "analytical boundary conditions" are proposed. (We intro- 
duce the term "analytical boundary conditions" as the boundary conditions that 
are required for the partial differential equations, so that the reader will 
not be confused with the term "numerical boundary conditions" that are required 
for the finite difference equations but not the differential equations.) The 
following is a summary of the conditions for well-posedness; refer to appen- 
dix B for a detailed derivation. In the following, we use k-j^j and tj_j as 
the ith row and jth column of the matrices and T“^, respectively, 

where and are defined as before (with frozen coefficient). The 

boundary is assumed to be at the left of the domain and the flow direction is 
from left to right. The gi’s and are given values. 

Subsonic inflow 0 < u < c : There are two positive eigenvalues of A. We 
require two analytical boundary conditions. The necessary and sufficient con- 
ditions for well-posedness are as follows. 

Conservative form: Impose any pair 



or impose 

'ki^p + k^2^ -f k^3e « gi(t) 

,^2lP ^ '^22™ *^23® “ S2^*^^ 

Non conservative form: Must impose p, that is, we have to impose 

(P /P /u 

I or I but not I 

lu Ip t P 

or impose 

'tjjP + + t^3P = g^CO 

.^2lP + + t^3P = g^Ct) 

Subsonic outflow (0 > u > -c) ; There is one positive eigenvalue. We 
require one analytical boundary condition. The necessary and sufficient con- 
ditions for well-posedness, for either the conservative or nonconservative 
form, are as follows: Impose any one of the variables or impose 
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kziP + + k^^e * g^Ct) 

for the conservative form, or impose 

t2iP + t22U + tjjP =* gaCt) 
for the nonconservative form. 

Supersonic inflow u > c : There are three positive eigenvalues. We 
require three analytical boundary conditions. 

Supersonic outflgy u < -c : There is no positive eigenvalue. Therefore 
we may not impose any analytical boundary condition. 


3. STABLE DIFFERENCE APPROXIMATIONS FOR HYPERBOLIC INITIAL 
BOUNDARY VALUE PROBLEMS (IBVP) IN A FINITE DOMAIN 


There are essentially three main considerations in studying approximate 
solutions to the initial boundary value problems (IBVP): (1) well-posedness 

of the original partial differential equations (PDE); (2) the method of con- 
structing extra boundary conditions required for the finite difference equa- 
tions (FDE) but not the PDE; and (3) the stability and accuracy of the FDE. 

In this section, we will review some of the well-known theory on stability 
analysis for IBVP, and list some of the commonly used stable schemes (stable 
for the combined interior and boundary schemes). The subject of accuracy will 
not be addressed here. The reader should refer to Gustafsson (refs 18, 19), 
Varah (ref. 20), Skdllermo (refs. 21, 22), and Sloan (ref. 23) for more detail. 
The major result for accuracy analysis is due to Gustafsson (refs. 18, 19), 
who proved that boundary schemes can be at most one order lower than the 
interior schemes, without loss of global accuracy. 

The treatment of difference approximations relating to Cauchy (initial 
value) problems of the hyperbolic type is quite well established. On the 
other hand, the treatment of mixed IBVP is considerably less well established. 
So far, the boundary conditions are quite difficult to construct and a poor 
choice can lead to inaccuracies and instabilities. The stability theory for 
difference approximations of the IBVP is really only complete for one space 
dimension, although this theory is essentially sufficient if the approximations 
are dissipative in the tangential directions (ref. 4) for multidimensional 
problems. For a one-space-dimension variable coefficient or quasi-linear sys- 
tem of hyperbolic equations with smooth solution (no shocks), the theor Is 
well established. Care is needed to avoid exponential growth due to imr per 
boundary extrapolation (refs. 9, 24). Recently Oliger (ref. 10) develop an 
easy way of constructing stable boundary schemes for the 1-D scalar prob 
For problems of higher dimension, little is known except for problems wii 
smooth boundaries, constant coefficients, and strictly hyperbolic cases. 

In the study of how boundary approximations affect the stability of gas- 
dynamics equations, rigorous stability analyses have only been applied to L-D 
and 2-D scalar equations with variable coefficients or quasi-linear property, 
or to systems of equations with constant coefficients. Boundary approximations 
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for problems with open boundaries and for viscous fluids at high Reynolds num- 
bers have not been studied sufficiently. Boundary approximations for factored 
or split ted implicit methods have not been analyzed. Crandall and Majda 
(ref. 25) have developed a complete treatment of the stability and convergence 
properties for scalar conservation laws in several space variables. Their 
method is a conservation-form, monotone difference approximation. Many inves- 
tigators have applied various boundary approximations to the nonconservative 
form of the nonlinear system and have compared the results with experimental 
data (see the Bibliography: Fluid Dynamics). Coughran (ref. 26) has devised 

a numerical method based on normal mode analysis (defined in appendix C) to 
study stable boundary schemes for the 1-D Euler equations. The following is a 
summary of the recent developments of currently available tools for stability 
analysis — concentrating on the more fundamental aspects of the subject, with 
a more detailed description of the theory for one space dimension. All of the 
initial data that we use throughout the report are assumed to be 
square-integrable . 


Fundamental Concepts 

In order to explain some of the difficulties, let us consider the differ- 
ential equation 


u^-u^=0 x>0,t>0 

u(x,0) = f(x) 

From the well-posedness of the problem, we know that no boundary conditions 
should be specified for x = 0, t > 0. If we want to solve equation (16), 
using some finite difference scheme, we need information about u at the 
"numerical boundary" x = 0, unless we use appropriate one-sided spatial dif- 
ferencing. For convenience, we will call the imposed boundary condition the 
"analytical boundary condition" and the extra boundary condition needed for 
the difference approximation the "numerical boundary condition," 

Let us say we want to solve the above problem using the leap-frog scheme 



n+i 

v 

J 


n-i 

V . 

J 




\ 


(17) 


where V4^ = v(jAx, nAt) denotes the numerical solution of u. We assume that 
At /Ax < 1 ; that is, equation (17) is a stable approximation for the Cauchy 
problem. We need an additional equation for v(0,t). Let us overspecify 
v(0,t) as 

v(0,t) = g(t) 

In general this overspecification will destroy the convergence. The only 
exception is the case in which v(0,t) = u(0,t), where u(0,t) denotes the 
solution at the boundary. But normally, we would not know about the exact 
solution. Kreiss and Lunqvist (ref. 27) and Gustafsson and Kreiss (ref, 17) 
have shown that "inexact" overspecification of boundary conditions leads to 
oscillatory solutions for this type of scheme (centered scheme, nondissipative) . 
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Therefore, one needs to be very careful when overspecifying boundary cor 
tions. The solution will look nicer if the approximation is dissipativ- 
because the oscillations will be damped. However, near the boundary th rors 
are quite serious. If one considers a system of equations, this error be 
propagated into the interior of the region by the ingoing characterist of 
the coupled variables, even when dissipative approximations are used, ,e 
stable way of handling v^^ is 


v 


n-i 

1 


Another is 


2v^‘ 


-1 


, n-2 


To illustrate another difficulty, let us consider 



u(x,0) « fi(x) 
v(x,0) - f 2 (x) 
u(-l,t) » g^Ct) 
u(l,t) = g^Ct) 

where the fi’s and g^'s are square- integrable. From the method of charac- 
teristics or normal-mode analysis, it can be shown that this problem is well- 
posed, In solving the equation numerically, we generally need special differ- 
ence equations to find v at both boundaries, even though analytically, the 
solution is uniquely determined for the PDE. Gottlieb and Turkel (ref. 24) 
have shown that if one uses the Lax-Wendroff finite difference method in the 
interior and quadratic spatial extrapolation for v at the boundary, then the 
resulting system is unstable, but Gustafsson et al, (ref, 9) have shown that 
the same extrapolation is stable in conjunction with the Lax-Wendroff method 
for scalar equations. In reference 28, Gottlieb et al. show that a straight- 
forward extension of the scalar results to a system may not work. However, by 
proper use of the characteristic variable at the boundaries, they demonstrate 
how the results of the scalar equation can be extended to a system. They show 
(in ref. 28) that by using quadratic spatial extrapolation for the appropriate 
characteristic variables, the revised method is stable. This is sometimes 
called the "characteristic stability theorem." 
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Overview and Development of Stability Theory 

For a one-space-diraension linear constant coefficient system, we can 
divide difference methods into two classes ~ those that belong to the "method 
of lines" approach and those that do not. The method of lines uses a finite- 
difference approximation in space and an ODE (ordinary differential equation) 
solver in time. For the method ^of lines approach, the stability of some of 
the popular schemes, like the central, forward, and backward spatial differ- 
encing schemes, coupled with simple boundary approximations was analyzed by 
Gary (ref. 29) (the matrix method), and by Dahlquist (ref. 30) (the positive 
real-function approach). In appendix D, we discuss the stability analysis of 
Gary and Dahlquist. They only showed the stability of the method for fixed 
that is, they did not show stability in the usual sense. In order to 
satisfy the definition of stability, these methods involve the additional 
analysis of infinite dimension matrices. 

For the approach that is not a method of lines approach, the simplest 
heuristic condition for stability was discussed by Trapp and Ramshaw (ref. 11). 
Their analysis used the interior as well as the boundary approximation to do a 
related Cauchy problem by the Fourier method (Von Neumann). An interior or 
boundary approximation is said to be Cauchy stable if it is stable for the 
related Cauchy problem (the related initial value problem , i. e . , the domain 
for 1-D is < X < °°) . They claimed that the minimum of the related Cauchy 

stability bound for the interior and the boundary can be used as the stability 
bound of the entire problem. But this heuristic approach does not provide 
sufficient conditions or proper hypotheses for stability of the IBVP. 

The most rigorous classical approach to the stability bound is the energy 
method. It is a powerful tool in dealing with certain particular equations 
or particular classes of equations (refs. 3, 6, 12). It can become rather 
complicated or tricky to apply, but it can deal effectively with boundary con- 
ditions and handles variable coefficients easily. However, it does not give 
necessary and sufficient conditions. 

A more unified approach to stability theory is due to Kreiss (ref. 8), 
and to Gustafsson et al. (ref. 9). It is sometimes called the normal-mode 
analysis. Strikwerda (ref. 31) has applied this theory for the method of lines 
approach. Godunov and Ryabenkii, whose work is discussed in reference 12, 
first gave necessary stability conditions for 1-D problems by considering 
modes of the form Uj^ (n - time step index, j - space mesh point 

index), where |y| <1 and j counts mesh points away from the boundary, Kreiss 
(ref. 8) and Gustafsson et al. (ref. 9) have greatly refined the approach, 
giving only mildly stricter conditions which are necessary and sufficient for 
stability. However, the analysis is more complex than that for the interior 
(i.e., the Cauchy problem). There are some important simple cases that have 
been studied in detail by this method, especially for dissipative approxima- 
tions. This theory is a posteriori in nature. Given a difference method, we 
can use this theory to determine whether the method is stable; but the stabil- 
ity criteria are often very difficult to verify. An example of how this theory 
applies to the first-order hyperbolic scalar equation with the simple well- 
known difference approximations can be found in appendix C. Recently, Oliger 
(ref. 10) gave sufficient stability conditions that are very easy to check. 
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A detailed discussion of the related theory is presented in appendix E. '^ese 
conditions can be used to guide us in the construction of stable methor r 
the initial boundary value problems. In order to make the development i 

understandable, we use (in appendix E) the case of a strictly hyperbol /s- 

tern with constant coefficients and coupled boundary conditions that ar ^11- 
posed. Then we discuss how we can arrive at the point at which it is c 

necessary to consider anything more complex than a single scalar equr n for 

each transformed variable. The stability analysis of this scalar eq^ Lon in 
a finite domain is equivalent to the analysis of two related quarter ane 
problems. We then proceed to discuss the way to construct stable schwines. 

The main assumption of the theory for constructing stable schemes is that the 
interior and boundary approximations are Cauchy stable and at least one of the 
approximations is dissipative. A point of caution — the sufficient condition 
does not guarantee sharp limits for conditionally stable methods. 


Some Stable Boundary Schemes (for Right Quarter Plane Problem, 

i. e. , x > 0) 


The following are some popular boundary schemes. 
Extrapolation; 


V 


v 


V 


V 


n+1 

o 


n+i 

o 


n+1 

o 


n +1 

o 


vn+i 

2^n+i _ ^n+i 

2vj^" - 


One-sided scheme: 



(18a) 

(18b) 

(18c) 

(18d) 


Box scheme: 



By using the normal-mode analysis (ref. 9), it can be shown that using the 
boundary schemes ((18a), (18b)), the one-sided scheme and the box scheme, 
together with the Lax-Wendroff or the Crank-Nicholson method, produces stable 
schemes for the right quarter-plane model problem: 


u^-Ujj»0 x>0,t>0 
u(x,0) - f(x) 


(19) 
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Extrapolation at the same time level ((18a), (18b)) (spatial extrapolation) is 
not a stable process for the leap^^frog scheme. For leap-frog types of schemes, 
we have to use (18c) and (18d), the one-sided scheme, or the box scheme. For 
predictor corrector schemes, like that of Richtmyer and Morton (ref, 12), or 
the MacCormick scheme (for linear constant coefficient, these two methods are 
identical) there are intermediate steps involved; Gottlieb and Turkel (ref, 24) 
have studied these schemes in detail. They have shown that spatial extrapola- 
tion ((18a), (18b)) and the one-sided boundary schemes are good choices. 

Now, we consider the class of interior schemes that evolves from linear 
multistep methods in ordinary differential equations (ref. 32), For model 
equation (19) with central spatial differencing, this class of schemes is of 
the form 

/u" , - u? \ 
p(E)uj" = -fltp(E){ 2Ax ^ / 

Here E is the shift operator defined by 

3 J 

and p and a are polynomials defined by 

P(E) = (1 + C)e" - (1 + 25)E + C 
a(E) = 6E^ + (1 - e + (|))E - 

The notation is consistent with that for linear multistep methods for 
ordinary differential equations and p(E) should not be confused with density. 
Some of the well-known methods (in time) belonging to this class are listed 
in table 1, 


TABLE 1.- PARTIAL LISTING OF LINEAR MULTISTEP 
METHODS 


f 

Method 

c 

e 


1 1 . Backward Euler 

0 

1 

0 

: 2. Two-step backward Euler 

-1/2 

1 

0 

; 3, Trapezoidal (Crank-Nicholson) 

0 

1/2 

0 

, 4. Backward differentiation 

1/2 

1 

0 

j 5 . Adams 

0 

3/4 

-1/4 

1 6. Lees 

-1/2 

1/3 

-1/3 

i 7. Two-step trapezoidal 

-1/2 

1/2 

-1/2 

i 8. A-contractive 

-1/6 

5/9 

-2/9 

; 9. Leap-frog 

-1/2 

0 

0 

1 10. Milne 

-1/2 

1/6 

-1/6 
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The first eight methods are unconditionally stable for the Cauchy problem and 
the remaining ones are conditionally stable. For the class of all two-step 
methods that are at least second-order time accurate, the parameters (9,C,4>) 
are related by (ref, 32) 


({j - C - 0 + 1/2 

The class of all third-order methods (in time) is obtained by imposing the 
additional condition 


^ = 28 - 5/6 

There is a unique fourth-order method, specified by 

9 3 -((> s -5/3 » 1/6 

which is called Milne* s method- Asume X = (At/2Ax) is chosen such that the 
method being discussed is stable for the related Cauchy problem. That is, 
equation (20) with j = 0, ±1, ±2, . . . , is stable. Gustafsson and Oliger 
(ref. 33) have proved the following results: 

I. If the boundary extrapolation ((18a) and (18b)) is used with the 
method (20) in table 1, then the resulting methods are all stable for the 
initial boundary value problem (19) except for the leap-frog and Milne methods. 

II. If the boundary extrapolation (18a) and (18b) is used with the 
method (20) in table 1, then the resulting methods are all stable for the 
initial boundary value problem (19) except for the two-step backward Euler, 
Lees, and two-step trapezoidal methods. 

All of the numerical schemes (interior ^ boundary) that we are going to 
study in the next section are mainly implicit schemes. For the model equa- 
tion (19), these schemes are unconditionally stable. One of the schemes is 
the backward Euler method in equations (20) and (18b). 


Stability Analysis of a Finite Domain 
Consider the scalar hyperbolic equation 

Ut 0<x<l,t>0 (21) 

with initial condition u(x,0) ■ f (x) . From the well-posedness of the problem 
(ref. 1), we have to specify analytically a boundary condition at the right 
boundary x ■ I, when c is negative, or at the left boundary x = 0, when 
c is positive- Hence, in addition to equation (21) and the initial data, ve 
specify boundary conditions 

u(l,t) = gi(t) if c < 0 
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or 


u(0,t) = go(t) if c > 0 

Let us assume c < 0. From a theorem of Gustafsson et al. (ref. 9), the 
stability of a difference approximation for the initial boundary value prob- 
lem (21) on 0 < X < I is equivalent to the stability of two related quarter- 
plane problems. The related right and left quarter-plane problems are defined 
as 


u 


t 


+ cu , 


u(x,0) 


0 

f(x) 


0<x<», t>0 
c < 0 


( 22 ) 


U^ + CU^ = 0 -°o<X<l,t>0' 

c < 0 

u(x,0) - f(x) 

u(l ,t) = g^(t) ^ 


(23) 


respectively. If only two- and three-point schemes are considered, then the 
stability analysis of the IBVP associated with (3.6) is transferred to the 
right quarter-plane problem (22). The stability of the left quarter-plane 
problem (3.8) reduces to the stability of a Cauchy problem. 


4. APPLICATIONS TO THE 1-D INVISCID EQUATIONS OF GAS DYNAMICS 


From the computational point of view, the unsteady inviscid gas-dynamics 
equations (Euler equations) in conservation law form have the following 
properties : 

a. They are a quasi-linear hyperbolic system, 

b. In general, the Jacobian of the flux vector consists of mixed posi- 
tive and negative eigenvalues (characteristic speed) . 

c. The flux vectors of the Euler equations are homogeneous functions of 
degree one in the dependent variables. 

d. The homogeneous properties provide a formal procedure for decomposing 
the flux vectors into subvectors, each of which depends on eigenvalues of the 
same sign (flux-vector splitting (ref. 16). Consequently, one-sided spatial 
difference operators can be used to construct a dissipative scheme. 

There are essentially two popular forms of the Euler equations being 
used in the computational fluid dynamics field; the conservative form (9) and 
the nonconservative form (11). Mathematically they are equivalent, but from 
the computational point of view they produce different solutions, if the same 
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numerical scheme is applied on the two forms* The study of well-posedness and 
stability of difference approximations is easier using the nonconservative 
form. Most of the existing theory and applications of the theory on the above 
studies use the nonconservative form. Recently, the development of computa- 
tional methods utilizing one-sided differencing gained popularity. Some of 
the one-sided differencing schemes are those of Godunov (ref. 34), Steger and 
Warming (ref. 16), Engquist and Osher (ref. 35), Roe (ref. 36), Carver 
(ref. 37), and Lax and Harten (private communication). At this time, there are 
no published results comparing the above one-sided differencing schemes, but 
some are more difficult to use than the others. The flux vector splitting 
method is useful for the application of a one-sided dissipative scheme on the 
conservative form, since the method is very simple to use and provides a proper 
way of handling the inflow-outflow boundary efficiently. For example, we can 
apply the one-sided difference operators on the split-flux subvectors over the 
interior and boundary points or, we can apply the one-sided difference opera- 
tors on the split-flux subvector over the boundary points only. 

We are going to discuss the stability of a few numerical schemes for the 
I-D Euler equations. Stability analysis is based on local linearization and 
solutions are assumed to be smooth near the boundaries. The various methods 
of handling the numerical boundary will be discussed briefly, but the method 
of linear extrapolation in the characteristic variables will be the main topic. 
Some numerical solutions of the quasi-l-D nozzle problem will be used to illus- 
trate the commonly discussed issues; for example, explicit versus implicit 
boundary schemes, unconditionally stable schemes, and underspecification or 
overspecification of boundary conditions. 


Flux-Vector Splitting 

As discussed earlier, the nonlinear flux vector F(U) is a homogenous 
function of degree one in U; that is F(aU) * aF(U). By application of 
Euler's theorem on homogenous functions, it follows that 


F can be split into two parts as (ref. 16) 

F - f'*’ + F" 

where f"*” corresponds to the subvector associated with the positive eigen- 
values X"*" of A, and F” corresponds to the negative eigenvalues X"*. 
Therefore 


with 


F - f'*' + F » (A'*' + A”)U 
Q = MT , Q~^ = T"^Q"^ 


and matrices 


A+ = QA+Q-^ , A" = QA"Q"^ 
and T“^ are defined in section 2. 


1 1 IT*— 1- 


(24) 


"I ii: 
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For 0 < u < c , we have 




2y 


2yu + c - u 
2 (y - l)u^ + (u + c)^ 


,(y - l)u^ + 


3 ^ (u + C)^ ^ (3 - Y)(u + c)c- 


2(y - 1) 


(25) 


F" = — 
2y 


u - c 


(u - c) 


(u - C)^ (3 - y)(u - c)c^ 

2 2(y - 1) 


(26) 


For u > c , we have 


F"^ = F , F" = 0 


The diagonal matrices A^, A are given by 


A^ = 


+ |ul 

2 

0 

0 


u + c + ! u + c I 
2 

0 


U - C + U - C ! 


A = 


2 

0 

0 


U C ~ I U + C I 
2 

0 


u - c - u - c 


Difference Approximations of the Inviscid Equations 
of Gas Dynamics 

By adopting the notation of Warming and Beam (ref. 38) and of Beam and 
Warming (ref, 32), the 1-D system of inviscid gas-dynamics equations can be 
approximated by a simple generalized three-level time differencing in the 
p (E) form as 

(l + uAt ^ a")d(E)u" = -At[a(E) - wp(E)](|J^ (27) 
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(Note: [I + inAtO/3x)A"]p(E)U'^ » p(E)U" + wAt{3 [A^P (E)U"] }/3x. ) The param- 

eter 0 ) *[0/(1+ C)] is determined by the particular time-differencing approxi- 
matlon used. Scheme (27) includes the following well-known implicit formulas 
(see sec. 3); 

5=0,9*y,(j)»»0 trapezoidal (Crank-Nicholson) 

C=0,9»l,(j)»0 backward Euler 

backward differentiation 

In (27), » A(U^) , (8F/3x)^ = [3F(U^)/3x], and is the solution at 

t - nAt with AT as the time step. 

There are two ways to utilize the flux-vector splitting; 

a. Apply one-sided approximations on the split-flux subvector throughout 
the entire computational domain of definition. (For example, use backward 
spatial differences for the ’’positive” subvector and forward differences for 
the "negative” subvector.) 


b. Apply the one-sided approximations on the split-flux subvector on the 
first and last interior points only. 

If we apply the flux-vector splitting on both the interior points and 
boundary points, system (27) can be expressed in the following form 


[l + .4t(iA-^"+iA-“)]o(E)u''.-att,(E) +(^)] (28) 

where A"*", A”, F"^, F" are defined as in equations (2A)-(26). One-sided 
first-order backward and forward-difference operators can be used for the 
spatial derivatives on the left-hand side of (28): 


3x 


3x 


[A+%(E)U'^] 


[A-"p(E)u"] 


A+%(E)Uj" - At^^p(E)u"_^ 


Ax 


Ax 


0(Ax) 


0(Ax) 


The spatial derivatives on the right of (28) can be approximated by the 
first- or second-order approximations. The second-order approximations are of 
the form; 


I I j 
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3F' 




9x 


3F+ - 4FT , + Ft 

— IZl Izl + 

2Ax ^ 


O(Ax^) 


3F~ 

9x 


-SF:"" + 4FT^ - FT^ 

— ^ + o(.x', 


The resulting algorithm (for the associated linearized Cauchy problem) is a 
dissipative, unconditionally stable, second— order algorithm, if we use two- 
step backward Euler time-differencing (6 - 1 , ? = 1/2 , (|> - 0) . The solution 
of (28) requires block tridiagonal inversion. We can introduce an approximate 
factorization of the left-hand side of (28) , and change system (28) to the 
product of two operators as 




+ wAt A"'')p(E)u" 


-At[a(E) 






(29) 


The solution of (29) only requires block bidiagonal inversion. The sta- 
bility of (29) is more difficult to analyze. We will only use form (28) for 
the quasi-l-D nozzle. 


Instead of using one-sided differencing throughout, we use system (27) 
without splitting A and F into two parts in the interior. The spatial 
derivative can be approximated by central differencing. For the first and 
last computational points, we can use the form (28), 

So far, stability analyses of variable coefficient or quasi-linear 
hyperbolic problems are only known for scalar equations or for systems with 
smooth coefficients and smooth solutions (ref. 2), For systems with nonsmooth 
coefficients or solutions, nonlinear instability can occur; for example, when 
an eigenvalue changes sign. One remedy is to use a dissipative scheme or add 
a dissipative term to the original differential equation. The one-sided 
spatial difference schemes comes with dissipation and frequently we have no 
control over it. The centered (spatial) schemes require ’’added" on dissipa- 
tion but allow different dissipative weight treatment in different regions of 
the solution. Both methods are quite popular in the computational fluid 
dynamics field. 


Stability Analysis 

As we have discussed before, theory for stability analysis of difference 
approximations for 1-D nonlinear hyperbolic equations has been established 
only for schemes that are dissipative or for problems with smooth solutions. 

The method of analysis depends on the "freezing method." If we freeze (dF/9U), 
then there is no distinction between the conservative and the nonconservative 
form. For each x = Xq , t = t^ we have a system of constant coefficient 
equations to analyze; that is. 
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( 30 ) 



with 



u(Xo.to) = Uq 

u(Xo,to) + cCXo.Cq) - Uq + Cq 

u(Xo.tQ) - c(Xo,to) - Ug - Cq 


Thus, at each time-step, the stability analysis consists of a pointvise 
examination of equations (31). For the higher order explicit methods, it is 
easier to use 01iger*s method (ref. 10) than the normal-mode analysis method 
to check for stability. On the surface, Oliger's sufficient condition con- 
sists of two parts (assuming the combined Interior and boundary schemes are 
stable for the model problem). 

a. Apply the interior difference scheme to (30) and do Cauchy stability 

shecks for all that are interior points. 

b. Apply the boundary difference schemes to (30) and do Cauchy stability 
checks for all Xq that are boundary points. 

If conditions (a) and (b) pass the stability tests at each point for 
every time step, what can we say about the stability of the original uncoupled 
nonlinear system? Stability is confirmed if at least one of the ap-^oxiraa- 
tions is dissipative (this is a sufficient condition; that is, an stable" 
boundary scheme for the related Cauchy problem does not imply that com- 
bined — interior plus boundary — scheme is not stable) and if the sc ions 

are smooth. In the actual case, the stability checks of part (a) anc ) 
involve scalar equations only. For popular numerical schemes, Cauchy ibi ' ty 
bounds are known. The major work is the testing of the values of 
i = 1, 2, 3 at each grid point and time step. This is trivial since 
are known. The method of normal mode analysis can follow the same appr ich, 
except in this case we have a necessary and sufficient condition. But Higher 
order methods are more difficult to verify. Often, we have to resort to 
numerical methods of solving a set of complicated resolvent equations 


ni ii: 
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(defined in appendix C) . For problems with shocks, there is no guarantee 
that stability of the ’’freezing’* family will imply stability of the original 
nonlinear problem. But, usually, it is quite promising if we use a dissipative 
scheme. 


The Numerical Boundary Conditions 

To simplify the discussion, let us assume that the spatial differencing 
we are going to use will be a first-order-one-sided or central-difference 
scheme, and denote the left and right boundary node index as 0 and J, Then 
the spatial differencing of (27) and (28) on the first and last computational 
points involves terms like 


E,AU." 

J J 

where Eq , Ej are some known matrices determined from the previous time step, 
and The AUq^, AUj^ are partially known from the analyti- 

cal boundary condition, with the exception of supersonic inflow, A few of the 
popular methods of obtaining the expression for the numerical boundary condi- 
tions are by 

a. Extrapolating in space or space-time (refs. 28, 39). 

b. Discretizing the Riemann invariant equations (the nonlinearized form 
of Che characteristic equations) or the characteristic equations (12) locally 
(refs. 40, 41, and J. Oliger (private communication)). 

c. Taking derivatives of the known condition in order to produce an 
extra boundary condition (refs. 19, 42, and M. Hyman (private communication)). 

d. Using nonreflecting boundary conditions (refs. 17, 43, 44). 

e. Overspecifying the boundary conditions. 

For implicit schemes, methods (a)-(d) above are quite complicated to 
implement into a computer code. Method (e) is of limited usefulness since it 
requires ^ priori knowledge of the exact solution to the difference equation 
at the boundary. Method (a) has the advantage of being the easiest to use; 
therefore, our study concentrates on method (a). But, as we know, extrapola- 
tion procedures suffer from the disadvantage of not modeling the differential 
equation (or not depending on the differential equations). However, if we use 
spatial linear extrapolation together with the two ways of utilizing the flux- 
vector splitting from the preceding subsection (Stability Analysis), the 
spatial differencing is already tailored to the direction of the characteris- 
tic curve locally. The extra unknowns that are required at the boundaries are 
due to the noniterative property of the scheme and the coupling of the physi- 
cal equations. Therefore, the numerical procedure for the extra unknowns at 
the boundaries should not be as crucial — spatial linear extrapolation appears 


31 



to be a good choice. As before we will use the term "numerical boundary con- 
ditions" as the extra boundary conditions that are required for the FDE but 
not for the PDE, or as the extra unknowns at the boundaries due to the non- 
iterative property of the scheme (local linearization). 

In implementing any methods (a)-(d), there are numerous and complicated 
details involved. Here, we will simply consider the spatial linear extrapola- 
tion in detail. The main point of this study is to show that the use of spa- 
tial linear extrapolation as boundary schemes for the implicit method (dis- 
cussed in the subsection "Difference Approximations of the Inviscid Equations 
of Gas Dynamics," sec. 4), is quite successful. Other comparisons of methods 
and application to different types of physical problems will be reported 
elsewhere (ref. 45). 


Spatial Linear Extrapolation for the Numerical Boundary Conditions 


For physical reasons, we sometimes prefer to specify boundary data in the 
primitive variables and compute in terms of the conservative variables in the 
interior. Thus the choice of variables for the analytical boundary conditions 
to be imposed and ntanerical (or extra) boundary conditions to be extrapolated 
for the conservative form (9) can be divided into the following four groups: 


Group 

.... . 

Variable 
(anal. B.C.) 

Variab le 
(nura. B.C.) 

I 

Conservative 

Conservative 

II 

Conservative 

Characteristic 

III 

Primitive 

Primitive 

IV 

Primitive 

Characteristic 


Under certain inflow-outflow combinations, not all of the above ways of 
imposing analytical boundary conditions are mathematically possible (or physi- 
cally desirable). If possible, group I is by far the simplest to implement 
with the rest appearing in increasing order of complexity. Group IV, on the 
other hand is more physically desirable and more theoretically sound (ref. 28), 
Groups II and IV reduce to the scalar model hyperbolic equations for the 
linearized equations of (9) and (11), respectively. We can have a whole class 
of stable schemes to choose from, as discussed in section 3. This is also 
true for group I in the supersonic inflow or supersonic outflow case. Now we 
turn to discuss group III. For the subsonic inflow case, it has been shown by 
Gustafsson and Oliger (ref, 33) that all the approximations (27) with param- 
eter values in table 1 are stable, with the following boundary conditions: 

Pq" given 


given 


„ n _ n _ n 

Pq - 2pi - P2 
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For the subsonic outflow case, Gustafsson and Oliger (ref. 33) also proved 
that all the approximations (27) with parameter values in table 1 (except for 
leap-frog and Milne) are stable, with the following ways of handling the 
boundary conditions: 


(i) 


(ii) 



given 

= 

= 

given 

= 


P 

P 


n 

2 

n 

2 


n 


n 


Here, we will describe the spatial linear extrapolations in the characteristic 
variables, that is, group II, Other groups can follow similar procedures. 

The relation between the conservative and characteristic variables is 


t t 


with U the vector of conservative variables, and W the vector of character- 
istic variables. The procedures for group II at inflow (left boundary) will be 

i) Make a first-order approximation: 


o o o 

il) Reorder into subvectors (U^)q" and (U^^)^'^ where (U^) " is 

the "analytical" boundary condition and (ull)^" is the "numerical" boundary 
condition. 

iii) Reorder into subvectors (W^)^" and (wll)^" where (W^)q" 

corresponds to the subvector associated with the positive eigenvalues of A 
and (W^I)q^ corresponds to the negative eigenvalues of A (for outflow right 
boundary, the signs of the eigenvalues are the reverse), 

iv) Reorder (T and partitioned it accordingly as 

Pi ?2 
P3 P4 
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Then we have 



Note that the delta foraulation (AU) is important for step (i) b luse of the 
nonlinear relation between conservative, primitive, and characteristic vari- 
ables. Now spatial linear extrapolation in the characteristic variables means 




This implies 

(PjAU^ P^AU^^)o" - 2(PjAU^ + P^AU^^)^” - (PjAU^ + P^AU^^)^" 

Since P^ should be nonsingular for a well-posed problem, we can rearrange 
terms and obtain 

(AU^^)g" - Rq(AU^)q" + Ri(AU)j" + Rj(AU)2“ (32) 

where Rx» R 2 known rectangular matrices which can be evaluated from 

the previous time step. (Note the mixture of dimensions in the equations.) 
Similarly, the outflow numerical boundaries can be expressed as 

(AU^^)j“ - S^(AU^)j" + S2 (AU)j_^ + S2(AU)"_^ (33) 

A similar formula can be derived if we impose the analytical boundary 
condition with the primitive variables (group IV) 

(AU^^)o” - Rq(A0^)q" + Ri(AU)j” + R^CAU)^" 

for the inflow boundary. By imposing primitive variables as analytical bound- 
ary conditions for the conservative system, group IV involves extra lineariza- 
tion and extra computations. 

If instead of using linear extrapolation for the numerical boundary con- 
ditions we discretize the characteristic equation and obtain an expression for 
(AU^^)q^, the counterpart of Rj[’s will be even more complicated than the 
Rl’s. 

There are two ways to alter the existing code by using the implicit 
boundary scheme! 

(a) Add correction matrices like (32) and (33) onto the first and last 
block rows of the block tridiagonal matrix. 
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(b) Use equations like (32) and (33) as extra equations in the block tri- 
diagonal matrix; that is, increase the dimension of the block tridiagonal 
matrix by dim(Uo^^) dim(Uj^^) — dim(Uo^^) Cleans the dimension of A 

word of caution, the final form of the matrix might not be in block tridiagonal 
form. 


Some Numerical Results 

The nozzles we consider are shown in figures 6 and 7 (refs. 39, 46). We 
use the unsteady gas-dynamics equations to obtain the steady-state solutions 
for various inflow-outflow conditions. The numerical spatial derivative 
approximations for the quasi-l-D nozzle problem are summarized as follows in 
table 2. The time differencing is the backward Euler method (high in stabil- 
ity). The trapezoidal formula, although yielding greater accuracy for small 
CFL numbers, results in instabilities for large CFL numbers. Additional time- 
differencing approximations and numerical boundary condition procedures will 
be considered in a future paper (ref, 45), 



A{X) = 1.398 + 0.347 ♦ TANH (0.8 X - 4) 

Figure 6,- Shubin nozzle, (ref. 46) for supersonic inflow, subsonic outflow 

study. 



A(x) — 1 + (A^j^ - 1) [(Xyi^ - X ) / Xji^] ^ ^ ^ ^TH 

A(x) = 1 + (Agx - 1) [(x - Xj,^) / (X^x - Xjh'J ^ > ^TH 


Figure 7.- Convergent-divergent nozzle (ref. 

outflow study. 


39) for subsonic inflow. 
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TABLE 2.- NUMERICAL SCHEMES 


Method 

Interior 

Boundary , numerical 

1 

Second-order one-sided^ 
(flux-vector splitting) 

Linear extrapolation 

2 

First-order one-sided 
(flux-vector splitting) 

Linear extrapolation 

3 

h Q 

Central •' + spectral norm 

(equivalent to Scheme 2) 
(ref. 16) 

Linear extrapolation 

4 

Central^^^ 

Linear extrapolation 

5 

Central + one-sided at 
first and last computa- 
tional points 

Linear extrapolation 


^Second-order for 9F"^/3x and 9F^/9x, but first-order 
for 3A'*’/3 x and 3A“/3x. 

^Fourth-order dissipation was added for the interior 
scheme. 

^Second-order dissipation was added at the boundary 
points. 


The numerical boundary conditions are treated either explicitly (E) , set 
to values at previous time step (replace n by n - I on the right-hand side 
of eqs, (32) and (33)), or implicitly (I), alteration of appropriate block 
tridiagonal matrix elements. 

The numerical scheme for each numerical experiment is defined by the 
temporal differencing (5,0,ct>), the spatial differencing (method 1, 2, 3. 4, 
or 5 of table 2), the variables chosen for the boundary conditions (groups I, 
II, III. or IV), and the temporal treatment of the boundary conditions 
(E or I). These choices obviously provide a large array of combinations which 
we must selectively sample. 

Typical steady-state solutions for three different flow conditions are 
shown in figures 8-10. Tables 3-8 present some of the results of numerical 
stability investigations. The calculations were made with a series of fixed 
CFL number and the numerical stability recorded. 

Although not extensive at this time, several general observations can be 

made: 


a. The results with boundary conditions I and II are very similar. 
Although the solutions are slightly different in the vicinity of the shock, 
the extrapolation of the conservative variables produces results that are 
comparable to those obtained when the characteristic variables are extrapo- 
lated (see tables 3-6) . 







Figure 8.- Density distribution; supersonic inflow, subsonic outflow, 

Shubin nozzle. 

b. For some schemes (see tables 7 and 8), the explicit and implicit 
treatment of the numerical boundary conditions produce similar numerical sta- 
bility bounds; that is, implicit treatment of numerical boundary conditions is 
not necessary for CFL > 1 (for some schemes). 

c. Overspecification of exact boundary conditions causes no problems. 
Figure II shows the supersonic inflow-outflow case. 

d. Methods 2 and 3 of table 2 behave almost identically. 

For the supersonic-subsonic problem, if we underspecify the boundary con- 
dition at the outflow, that is, without specifying anything, the solution 
diverges. Moreover, updating the boundary points via the delta form (32) 
and (33), and then obtaining 


AU " + U ” 
o o 
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VELOCITY 



Figure 9.- Velocity distribution; subsonic inflow, subsonic outflow, 
convergent -divergent nozzle, area ratio 2:1.16. 
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DENSITY 



X 


Figure 10.- Density distribution: subsonic inflow, subsonic outflow, 
convergent-divergent nozzle, area ratio 2.5: 1.5. 


TABLE 3.- NUMERICAL STABILITY CHART: 

BOUNDARY SCHEME I, SHUBIN NOZZLE 
(Boundary conditions: inflow = c,Tn,e; outflow = p) 

! I I I 

I j Method ! Method Method 

! i 2 ; 4 ^ 5 


(I) Yes : (I) Yes | (I) Yes ' 
(E) Yes j (E) No : (E) Yes ? 


I 

(I) 

Yes 

(I) 

Yes 

(1) 

Yes 


(E) 


Yes 



(E) 

Yes 

10 

(I) 

Yes 

(I) 

Yes 

(I) 

Yes 

(E) 

Yes 

— 

— 

(E) 

Yes 

! 

100 j 

' (I) 

No 

(I) 

i 

No ) 

(I) 

No 

(E) 

No 

— - 

! 

1 

(E) 

No 


Notes: Supersonic inflow; subsonic outflow. 

I = implicit numerical boundary condition 
E = explicit numerical boundary condition 
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TABLE 4.- NUMERICAL STABILITY CHART: 
BOUNDARY SCHEME I, SHUBIN NOZZLE 
(Boundary conditions: inflow = p,m,e; outflow » 


CFL 

Method 

2 

Method 

4 

Method 

5 

1 

(I) Yes 
(E) Yes 

(I) Yes 
(E) Yes 

(I) Yes 
(E) Yes 

5 

(I) Yes 
(E) Yes 

(I) Yes 
(E) No 

(I) Yes 
(E) Yes 

10 

(I) Yes 
(E) Yes 

(I) Yes 

(I) Yes 
(E) Yes 

100 

(I) No 
(E) No 

(I) No 

(I) No 
(E) No 


Notes: Supersonic inflow; subsonic outflow, 

I = implicit numerical boundary condition 
E * explicit numerical boundary condition 


TABLE 5.- NUMERICAL STABILITY CHART: 

BOUNDARY SCHEME II, SHUBIN NOZZLE 
(Boundary conditions: inflow « p,m,e; outflow = p) 


CFL 

Method 

2 

Method 

4 

Method 

5 

1 

(I) Yes 
(E) Yes 

(I) Yes 
(E) No 

(I) Yes 
(E) Yes 

5 

(I) Yes 
(E) Yes 

(I) Yes 

(I) Yes 
(E) Yes 

10 

(I) Yes 
(E) Yes 

(I) Yes 

(I) Yes 
(E) Yes 

100 

(I) No 
(E) No 

(I) No 

(I) No 
(E) No 


Notes: Supersonic inflow; subsonic outflow. 

I « implicit numerical boundary condition 
E = explicit numerical boundary condition 
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TABLE 6,- NUMERICAL STABILITY CHART: 

BOUNDARY SCHEME II, SHUBIN NOZZLE 
(Boundary conditions: inflow = p,m,e; outflow - m) 


CFL 

Method 

2 

Method 

4 

Method 

5 

1 

(I) Yes 
(E) Yes 

(I) Yes 
(E) No 

(I) Yes 
(E) Yes 

5 

(I) Yes 
1 (E) Yes 

(I) Yes 

(I) Yes 
(E) Yes 

10 

1 (I) Yes 
1 (E) Yes 

I 

(I) Yes 

(I) Yes 
(E) Yes 

100 

1 (E) Yes 
(E) No 

(I) No 

(I) No 
(E) No 


Notes: Supersonic inflow; subsonic outflow. 

I = implicit numerical boundary condition 
E = explicit numerical boundary condition 


TABLE 7.- NUMERICAL STABILITY CHART: 

BOUNDARY SCHEME IV, SHUBIN NOZZLE 
(Boundary conditions: inflow = p,u,p; outflow = p) 


CFL 

Method 

1 

Method 

2 

Method 

4 

Method 

5 

I 

(I) Yes 
(E) Yes 

(I) Yes 
(E) Yes 

(I) Yes 
(E) No 

(I) Yes 
(E) Yes 

5 

(I) No 
(E) No 

(I) Yes 
(E) Yes 

(I) Yes 
(E) No ; 
1 

(I) Yes 
(E) Yes 

10 

— 

(I) Yes 
(E) Yes 

(I) Yes ' 
(E) No 

(I) Yes 
(E) Yes 

100 

— 

(I) No 
(E) No 

(I) No 

(I) No 
(E) No 


Notes: Supersonic inflow; subsonic outflow, 

I = implicit numerical boundary condition 
E = explicit numerical boundary condition 
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DENSITY 


TABLE 8.- NUMERICAL STABILITY CHART: BOUNDARY SCHEME IV, 

CONVERGENT-DIVERGENT NOZZLE 
(Boundary condition: inflow = p, p; outflow » p) 



Notes: Subsonic inflow; subsonic outflow; 

area ratio 2:1.16. No shock- 
I =» implicit numerical boundary condition; 
E * explicit numerical boundary condition. 


i 

I 



0 2 4 6 8 10 

X 


Figure 11.- Density distribution; overspecify at outflow, 

Shubin nozzle. 
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where the R^'s and S^'s are the same as in equations (32) and (33), produces 
a solution that is not as smooth near the boundary. 

The smoothing parameter for the fourth-order dissipation term for 
methods (3)- (5) of table 3 are 0.5, No study has been made for varying the 
smoothing parameters for different solution behavior zones. 


CONCLUSIONS 


A comprehensive overview of the state of the art of well-posedness and 
stability analysis of FDE for IBVP of the hyperbolic type was presented. The 
^freezing'* theory was used as a guide to construct boundary schemes for the 
1-D inviscid gas-dynamics equations. The use of primitive variables as the 
analytical boundary conditions for the conservative form of the 1-D inviscid 
gas-dynamics equations was formulated and then applied to the quasi-l-D nozzle 
problem. 

Spatial linear extrapolation as a boundary scheme can produce reasonable 
steady-state solutions. It is scheme-independent, and thus provides a compact 
form for computer code implementation. Added dissipation terms, the linear- 
ization of the (3F/dU) matrix and ways of updating the boundary points can 
affect the stability and accuracy of the solution. Future work in this area 
is needed. 
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APPENDIX A 


DEFINITIONS AND EXAMPLES OF WELL-POSED HYPERBOLIC 
DIFFERENTIAL EQUATIONS IN '’Lj NORM” 


Well-Posedness of Cauchy Problem 

Definition A.l: The L 2 norm of a vector function u(x) with -00 < x < <» 

is defined as 


lu(x)|| = [jr; u* (x)u(x)dxj 


1/2 


where u* is the transpose and complex conjugate of u. 


Consider the Cauchy problem 

+ Au^ + Bu=*0 -«<x<‘»,t>0 

u(x,0) = f(x) 

where A and B are N x N constant matrices » and u and f are vectors with 
dimension N. 


(Al) 


Definition A. 2: For all initial values f(x) with l|f(x)j] < the 

Cauchy problem (Al) is well-posed if there are constants k, a (independent of 
f(x)) such that for ail solutions and all t, there exists an estimate. 


where 


u(x,t)|| < K e^^ ||u(x,0)|| 


u(x,t)|l 



u*(x, t)u(x,t)dx 


1/2 


(A2) 


Example: Consider the first-order scalar equation 

u^(x,t) - au(x,t) »0 -«»<x<«», t>0^ 

u(x,0) - f(x) j 

with ||u(x,0)|| < ® and known real constant a. The solution of (A3) is 

u(x,t) = f(x)e^^ 


(A3) 
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Since 


||u(x,t)|| = I |lu(x,0)|| 

the solution of (A3) satisfies (A2) . Therefore (A3) is well-posed. 

The definition of hyperbolic and strictly hyperbolic systems is as fol- 
lows: The system (Al) is hyperbolic if A is diagonalizable and with real 

eigenvalues. It is strictly hyperbolic if all the eigenvalues are real and 
distinct . 

There is a simple equivalent algebraic condition for definition (A, 2) to 
hold. This condition (ref. 2) is found by Fourier transforming equation (Al) 
in X and studying the norm of the Fourier transformed variable. Through 
this method, it can be shown that a hyperbolic system (Al) with all 
;iu(x,0)|| < and B = 0 are well posed. 

Let us define P(i'^) = -i-A with ua real. Then the algebraic condition 
is: The Cauchy problem for (Al) is well-posed if and only if there are con- 
stants K and a such that 

i P(i^) 1 1 _ at 

max 

Example: For the scalar hyperbolic equation 

u H- cu =0 

t X 

u(x.O) = f(x) 

with c real and ||u(x,0)i| < We have 

P(iuj) = -i^c 

* i c 

thus max e ^ <1. If we take K = 1, a * 0, then the algebraic condition 

is satisfied. That is |ju(x,t)|| =[|u(x,0)][. For the hyperbolic system (Al), 
the well-posed algebraic condition is immediately satisfied since there is a 
unitary matrix T s.t. 
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Well-Posedness of Initial Boundary Value Problem (IBVP) 


Consider the IBVP of the strictly hyperbolic system in the quarter plane 

(0 < X < 


+ Au^ » F(x,t) 
u(x,0) - f(x) 


0<x<«, t>0 


u^(0,t) = Su^^(0,t) + g(t) 

where A is an N x N diagonal constant coefficient matrix with 


(A4a) 

(A4b) 


(A4c) 


A" 0 


0 -A* 


j • 


u^) 


where Xj, j = I, , . N are real and distinct, S is an Jlx(N - Z) matrix, 
and f(x; is smooth. (It is no restriction to assume that A is in diagonal 
form because the system is strictly hyperbolic and can always be written in 
this form after a suitable transformation.) For simplicity, we will consider 
the homogeneous initial data u(x,0) =* f(x) * 0. The asstmption of homogeneous 
initial data is no restriction since we can always subtract that solution of 
the nonhomogeneous Cauchy (initial value) problem and obtain exactly this 
situation. 

Definition A. 3: We will say that the quarter-plane problem (A4) with 

homogeneous initial data is well-posed if the estimate 

||u(x,t)||^ dt < ^ 
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holds with a constant Independent of g and F, but perhaps depending on 

T. In here, the Lz norm of u(x,t) is defined as 


l!u(x,t)i| J u*(x,t)u(x,t)dx 

We can extend this definition for higher dimension systems of equations 
with a slight modification. Consider the IBVP of a two-dimensional strictly 
hyperbolic system (see sec. 2 for definition) in the quarter plane 
(0 < X < *» -« < y < w) 

+ Au^ + BUy » F(x,y,t) 0 < x < * 1 

u(x,y,0) » 0 -« < y < > (A6) 

u^(0,y,t) - Su^^(0,y,t) + g(y,t) t > 0 j 

where A, S, and are defined as before (with A replaced by -jA + - 2 ®> 

where and are real and ■ 1). 

Definition A. 4; We will say that the quarter-plane problem (A6) is 
well-posed if the estimate 


f i;u(0,y,t)j|^ dt + f ;lu(x,y,t),;' dc 

1 /g(y.t)'J dt + J* dt^ 

holds. Here depends on t but not F and g. Where the Li norm 

are defined as 


liu(o.y.t)i!y 



u*(0 ,y, t)u(0,y ,t"^dy 


and 


||u(x,y,t)l|^ ^ “ r f ^*(x,y,t)u(x,y ,t)dx dy 

with similar definitions for j!g(y»t)i!y and l!F(x,y ,t)|ix,y* 

For the one-dimensional systems, we can get the same ccnditicns as in 
definition (A3) by using the method of characteristics. Tliis is not the case 
for higher dimensional systems (refs. 2, 6). The application of the method of 
characteristics is discussed in detail in section 2, Here wo will state the 
necessary and sufficient algebraic conditions for definition (A3)» This is a 
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simplified version of the main theorem of Kreiss (ref. 1). The theory of 
Kreiss (ref. 1) treats problems inany number of space dimensions. Interested 
readers should refer to reference 1 for extension to more space dimensions. 

For the two—dimensiohal and three-dimensional quasilinear systems of inviscid 
gas dynamics equations, please refer to Oliger and SunstrOm (ref. 6) and Oliger 
(ref. 7). 

The so-called "normal-mode analysis" algebraic conditions for defini- 
tion (A3) will be stated after the following brief preliminary background. Let 
us Laplace transform (A4a) and (A4c) with respect to t and denote s =* n + i? 
as the variable dual to t. We obtain 


su + Au^ = F 

for 

X > 0 

u = Su + g 

for 

X = 0 


The symbol (") is the Laplace transformation of the variable ( ). 

Associated with (A7) is the following eigenvalue problem. A square- 
integrable function <)>(x) for 0 < x_< °° is an eigenfunction of (A7) corre- 
sponding to an eigenvalue s if ({> is a solution of the problem 


Sif + A(j)^ 

- 0 

for 0 < X < 

(A8) 


= 

for X = 0 

(A9) 


We do not want s with ii = Re(s) >0 to be an eigenvalue of (A8) and (A9) . 

If this happens, is not in L 2 (<^ is not in L 2 means (j) is not square 
integrable). Therefore, we have to decide whether s with Re(s) >0 is an 
eigenvalue or not. Equation (A8) is an ordinary differential equation whose 
general solution in L 2 for Re(s) > 0 can be expressed as a linear combina- 
tion of . a linearly independent normalized eigensolutions (see Kreiss, ref, 1, 
for details). That is, the general solutions in L 2 depend on Z free 
parameters a » ...» Introducing the solution into (A9) , we get 

a linear system of equations: 

R(s)a = 0 , R a matrix function of s (AlO) 

and s is an eigenvalue if and only if 

Det[R(s)] =0 ^ (All) 

Kreiss has shown that Det[R(s)] is a continuous function of s for 
Re(s) > 0 and he defines s = i^ to be a generalized eigenvalue if 
Det[R(iC)] = 0. 

Now we can state the necessary and sufficient conditions for the estimate 
of type (A5) to hold. 


Theorem A. 1: The IBVP for (A4) is well-posed in L 2 if and only if the 
eigenvalue problem (A8) and (A9) has no eigenvalue or generalized eigenvalue 
for Re(s) > 0. 
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Next, we want to give an example to show that for a one-dimensional sys- 
tem, using the method of characteristics is equivalent to Theorem A, I, Here we 
assume that readers are either _ familiar with the subject of the method of 
characteristics or will consult references 2-4 (or sec. 2) for details. Con- 
sider the following quarter-plane problem for the wave equation 

w^t = ^xx 0 < X < °°, t > 0 (AI2) 


with initial conditions 


w(x,0) = f(x) 
w^(x,0) = g(x) 


and boundary conditions 


w(0,t) = w^(0,t) = 0 

We can recast the problem into a system of first-order hyperbolic form by 
letting V = w^ , u = w^^, = v - u, and Z2 = v + u. Then (A12) becomes 


^1 o\/z. 



= 0 


(A13) 


1/\Z 


with initial conditions 


'g - fx' 
<g + fx/ 


at t = 0 


(A14) 


and boundary conditions 




at X = 0 


(AI5) 


From the method of characteristics, we can see that the initial condition (A14) 
together with the boundary condition (A15) determine the solution of (A13) 
uniquely. 

Now we turn to Theorem A.l. The general solution of the associated eigen- 
value problem for (A13) in L 2 with Re(s) >0 is 


4) = Oi e*®"" 6i 

where e^ is the normalized eigenvector of 


(A16) 


corresponding to Re(s) > 0, The normalized eigenvector e^^ is found to be 



Introducing (A16) into (A15) , we get 




-sx 


x=»o 


= 0 


^50 


Therefore, there are no non-trivial solutions in Lz for Re(s) > 0 and the 
problem is well-posed. 
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APPENDIX B 


CONDITIONS ON WELL-POSEDNESS OF THE INVISCID EQUATIONS 

OF GAS DYNAMICS 


Freeze the coefficients of the Jacobian matrix A and rewrite equa- 
tions (10) and (11): 




|^+Ayi = 0, U = M"^U 

dt 9x t t 


and the characteristic equation 


+ A = 0 

dt dX 


with 


W = T~^U 


(Bl) 

(B2) 


(B3) 


(B4) 


or 


W = 


(B5) 


where U, U, W are the conservative, primitive and the characteristic vari- 
ables of (Bl) , (B2) , and (B3) , respectively. The matrices and T“^ are 
defined as 
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and Uq, Cq, Pq are the values of u, c, p due to freezing of coeffic :s. 

We want to discuss two sets of well-posed boundary conditions: ( 1 ) imr i 

boundary conditions (case 1) that are in the form of the individual na 1 

physical variables alone, that is, the primitive or the conservative ‘ if; 

and (2) imposed boundary conditions (case 2) that are in the form of inear 

combination of the physical variables, that is, a^^p + a2m + a^e = g(t or 

b^p + b2U + b^p “ g(t) where the b^’s, g(t) and g(t) are the .ven 

known quantities, not p, m, e, u, or p. 


Case 1 

Assume that we want to impose the analytical boundary condition in terms 
of the conservative variables. (We use the term "analytical boundary condi- 
tions" as the boundary conditions that are required for the partial differen- 
tial equation.) The boundary is assumed to be at the left of domain and the 
flow direction is from left to right. Thus the number of positive eigenvalues 
(equal to the number of analytical boundary conditions) and negative eigen- 
values is known. The procedure to check for well-posedness consists of two 
steps. First we reorder (B 5 ) as 



where W^ and W^^ are the characteristic variables corresponding to positive 
and negative eigenvalues of A, respectively. And u^ corresponds to the 
proposed analytical boundary condition variables and u^^ represents the rest 
of the variables. Second, we have to check whether exists or is 

emoty. Thus the necessary and sufficient condition for well-posedness is 
Q4 exists or is empty. 

Similarly, if we want to impose the primitive variables as analytical 
boundaries, we can reorder (B 4 ) as 

/w^\ /Qi Q 

\W^V \Q 3 Q 

where well-posedness here means exists or is empty. 

Therefore, under a type of inflow-outflow condition, once we have decided 
on a set of analytical boundary conditions, the way to check for well-posedness 
of (Bl) or (B 2 ) is to see if the determinant of (Q4) or (Q4) is equal to zero 
or not. The following are the determinants of Q4 and (if it is not empty) 
for various choices of inflow, outflow conditions. Again, we want to empha- 
size that the boundary is assumed to be at the left of the domain. Therefore, 
we only^need to investigate the determinant of Q4 (or 0 ^^), The form of a 
Q4 (or Q ) depends on how we order the variables in w^^ and u^^ (or u^^), 
which differ by a change of rows and columns or both; but the absolute value 
of the determinants are the same. 
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Pure supersonic; (u > c for inflow and ~c < u < 0 for outflow): 


Primitive variables: inflow — There are three positive eigenvalues. We 

require three analytical boundary conditions. Thus is empty. 

Primitive variables: outflow — There are no positive eigenvalues. We 

do not have to impose any analytical boundary condition. Thus^ Qi^ « T"^ and 
Det(Qj^) =Det(T“^). Note that Det(Qi^) means determinant of 

Conservative variables: The situation is the same as in the case of 

primitive variables. Therefore, the well-posedness conditions are to impose 
all three variables for the supersonic inflow case and none for the supersonic 
outflow case. 


Subsonic outflow: (-c < u < 0) : There is one positive eigenvalue. We 

require one analytical boundary condition. Therefore, we can propose the fol- 
lowing three choices. 


Primitive variables: analytical boundary condition — p 


Q4 




Det(Q^) 




Primitive variables: analytical boundary condition — u 


Q4 




1 

/2p _c_ 
’^o o 



Det(Q^) = — ^ 


Primitive variables: analytical boundary condition — p 

1 



/2 /2p, 


Det(Q^) 


-1 




Thus, imposing any one of the variables p, u, or p will result in a well- 
posed condition. 
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Conservative variables: analytical boundary condition — e 

(Y - l)Uo 


1 


(Y - 


2^^ 


Uq (y - 1)Uq^ 


-1 


(1 - Y)u, 


/2p^ y^p„c 


o • --'O'O 
2 /■-. 2 


Det (QJ 


^ (y - l)\io (y - Duj, 


■^»0 




Conservative variables: analytical boundary condition — m 




(Y - Duq^ 
2Co^ 

(y - l)u^ 
+ 

2^Po<=o 


- ( y -V 


(Y - 1) 
'^Po'^o 


Det(QJ 


(Y - 1)(Cq + Up) 


Conservative variables: analytical boundary condition — p 


(y - Du, 


-(Y - D 


Qu = 


-1 


(1 - Y)u, 

/Iq c 
o o 


Co 


(y - D 
’^PoS- 


Det(QJ 


-(Y - 1) 
^Po'^o' 


Again, for well-posedness , imposing any one of the variables p, m, or e 
result in a well-posed condition. 

Subsonic inflow: (Q < u < c) : 

Primitive variables: analytical boundary conditions — u, p 

» 0 « Det(Q^) . , , not well-posed 
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Primitive variables; analytical boundary conditions — p, p 

= — = Det(Q^) 

Primitive variables: analytical boundary conditions — P» u 

;; 1 


*^Po<^o 


Det(QJ 


In this case imposing u and p will produce an ill-posed problem . 
Conservative variables: analytical boundary conditions — m, e 


(Y - l)u. 




Det{Q^) 


Conservative variables; analytical boundary conditions — p, e 


Qu = 


-1 

‘'^p. 


(1 - Y)u, 

’^Po'^o 


Det(Qj 


Conservative variables: analytical boundary conditions — p, m 

^ = Det(Q^) 

•^Po^o 

In this case, imposing any pair (p,m), (p,e) or (m,e) will result in a 
well-posed condition. From the above examination, the only analytical boundary 
condition set that produces an ill-posed problem is (u,p) for the subsonic 
inflow case. 


Case 2 

In this case, we only can impose the characteristic variables correspond- 
ing to the positive eigenvalues of A (or A) in order to obtain a well-posed 
condition. 

For supersonic inflow, we can specify all three characteristic variables 
, w^ , and w^ , that is. 
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yT '^PqS 


P = l,(t) 


1,1 - , X 

u + p = g,(t) 


^ / 2 p 


o o 


for the primitive-nonconservative form (11) or conservative form (9), and 


1 - 


(Y - 
2 c„^ 


Uo (Y - Du^n 

I — ... 


" ^ (1 - y)u, 

y^Po 2/2PqCq _ 

P + 

.>^Po ’^Po^'o 

■ % 

- 


r -1 a - Y)u, 

•^Po 2/2 PqCq _ 

P + 

j- 

_/2p^ '^PqCo 


- 1 -* 

1 


(y - 1) 
.'^Po‘'o 


e » gi(t) 


e = g,(t) 


e = g,(t) 


for the conservative form (eq. ( 9 )) where g^'s and g^'s are the values which 
are supposed to be specified. 

For subsonic inflow, we only can specify w^ and Wa » that is. 


p r p ” g,(t) 


/2 


u 


+ —r — p = 82 (t) 
’^Po^'o 


for system (11) or (9) and 



for system ( 9 ), Again the s and g^*s are the values which are supposed to 
be specified. 
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APPENDIX C 


EXAMPLES OF THE THEORY OF GUSTAFSSON, KREISS, AND SUNDSTROM 
(NORMAL MODE ANALYSIS) 


Here we briefly review the stability theory of Gustafsson, Kreiss, and 
Sundstrom for the initial boundary value problem of the hyperbolic type for 
the leapfrog method. Please refer to their original paper (ref. 9) for more 
details . 


Consider the following equation 

^+c|^=0 0<x<l, t>o) (Cl) 

at 3x " “ I 

u(x,0) = f(x) / 

In addition to equation (Cl), we specify boundary conditions 

u(0,t) = g^(t) if c > 0 

u(l,t) = g^(t) if c < 0 

But, numerically, one needs boundary conditions at both x = 0 and x = 1. 
Therefore, a separate procedure is used to determine the numerical boundary 
conditions . 

Let us solve (Cl) by the leapfrog scheme with At as the time-step, and 
Ax as the mesh spacing. We will use the notation 

« v(jAx,t) - u(jAx,t) 
t = nAt 

Assume for the moment c=-l,0<x<°°, and approximate (Cl) by 


n+i 

V . 

J 


n-i 

V . 

J 


^ (v" 
Ax ^ j+i 




Vj° = f(jAx) 


and the numerical boundary condition at x = 0 by 


V 


n 

1 


(C2a) 

(C2b) 


(C3) 


The Gustafsson et al. stability theory for this case seeks a general solution 
of (C2a) and (C3) of the form 


57 



for appropriate complex scalars z and k. This substitution is made in both 
the difference scheme in the interior and on the boundary. The basic scheme 
(C2a) is assumed stable for the Cauchy problem, that is, (At/Ax) = A < 1 for 
the interior points. 

By substituting = z’^j (x) into (C2a) and (C3) , we obtain 


z^ - 1 » 


'j+i Az 


- 0 j - 1, 2, . . . 


^o “ 


(C4a) 

(C4b) 


Equation (C4) is defined as the resolvent equation . 

Letting the solution of (C2) be z^k'^ , we obtain the characteristic equa- 
tion for (C2) as 


K 


2 



1 = 0 


(C5) 


A necessary and sufficient condition for stability of the IBVP is that 
(C4) have no nontrivial bounded solutions 


j = o J 

with !z| > 1. An eigenvalue to the associated equation (C4) is defined as a 
nontrivial solution to (C4) with bounded vj = <^((<1 < 1) and |z| > 1. A 
generalized eigenvalue to the associated (C4) is defined as a nontrivial solu- 
tion to_(C4) with Vj = , and [k| =1 and jz| = 1, such that all solu- 

tions z,K of (C5) with |z| > 1, and sufficiently close to z and k, have 
|k| < 1. The equivalent necessary and sufficient condition for stability is 
that the associated (C4) have no nontrivial eigenvalues or generalized 
eigenvalues. 

The stability analysis consists of the following four stages: 

1. The order d (d = 2 in this case) of the resolvent (difference) 
equation (C4a) determines the general solution of vj (x) — a linear combina- 
tion of d solutions 

Vj = + 0^X2^ + . . . + 

where the x^^'s are the roots of (C5). 
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2. The root structure of (C5) determines the type of solution for vj (x) . 
In this case the roots of (C5) have the following properties (see ref, 9 for 


detail). If 

1 z 1 > 1 , then 1 < 1 1 < I 

and [<21 

> 1. 

This is an immediate 

consequence of 

the Cauchy stable scheme 

of (C2a) 

. If 

z 

= ei“, then 


!<i! > kal > 1 

for 1 

sin 0 1 

> 

X 


1 < 1 = j K 2 1 ~ 1 

for 1 

sin 9 1 

< 

X 


K ^ “ -1 , <2 * ^ 

for 

e 

= 

0 


< ^ = 1 y <2*-! 

for 

0 


r 


K, = <2 = ±i 

for 

sin 8 

- 

±X 


3. The assumption that the interior scheme (C2a) is Cauchy stable helps 
delete the unbounded solutions of vj — all solutions with \<i\ ^ • The 

theory says that the general bounded solution of (C4a) is then 



From lemma (5,1) of reference 9, only one root of the quadratic (C5) has 
modulus less than one. When |z| = 1, one or both of the roots of (C5) may 
have modulus one. If this is the case, the for the general bounded 

solution of (C4a) is defined by continuity to be that root which is the limit 
of the root k(z), |<(z)j < 1 for |z| > 1, as |z| 1. Thus 


4. After substitution of vj = kJ in (C4b) , if there exists a nontrivial 
bounded solution for jz| > 1, the difference schemes (C2a) and (C3) will be 
unstable. In this case 

(k - 1) - 0 (C6) 

Therefore, when k * 1, (C5) and (C6) have a nontrivial solution with z = -1. 
From item (2) above we know that this is a generalized eigenvalue and thus 
stability is violated. 

In many instances, the root structure of the characteristic equation (C5) 
is difficult to analyze. Another way of testing for generalized eigenvalues 
is as follows: 


With z = -1, we want to find out whether k « 1 is or We 

therefore make a perturbation calculation, and study (C5) in the neighborhood 
of z = -1. Let z=-l-5, 6>0 and < = 1 + e with 5,e small. From 
(C5) we get 
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K = 


- 1 ^ ■ 

l/Zz^ - ] 


Xz 

V\ Xz 

7 


2 


(1 + 6)2 - 1 ^ 1/1 

f(l + 6)2-11 

2 

-X(l + 6) - V 

IT -X(l + 6) . 

+ 4 


Since 5 > 0, at least one of the i » 1,2 is negative, and < = 1 + e 

is not K . Therefore z =» -1 is a generalized eigenvalue and thus 

stability is violated. 


Now, consider using 

^ n+i n ' 

* V, 


(C7) 


instead of (C3). The equivalent of equation (C6) becomes 

(z - ic^)Cj = 0 


For jic^l < 1 and !z| > 1 

|z - kJ >0 

Thus, (C2a) and (C7) constitute a stable difference method for the right half- 
plane problem. 


Stability of some other explicit and implicit schemes^ using the above 
approach, can be found in Oliger (ref. 15), Gottlieb and Turkel (ref. 24), 

Sloan (ref. 23), and SkCllermo (ref. 21). For multistep schemes, the stability 
criteria of this method are often very difficult to verify. Here, we are 
going to discuss an unconditionally stable scheme in which we use it for the 
quasi— 1-D nozzle problem. Let us solve (Cl) by backward Euler in time and 

difference in space . The numerical boundary condition at x » 0 is 
by linear extrapolation 


n+i n 

V . - V. 

J j 


n+i 

^o 


v/^n+i n+i\ 

- "j-i) ■ 


n+1 

r 

1 


n+1 

^2 


The characteristic equation of (C8a) 

k(z - 1) = Xz(k^ - 1) 
and the boundary scheme (C8b) satisfies 


At 

2Ax 


(< - 1)- 


1 


(C8a) 

(C8b) 


(C9) 


(CIO) 
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The Cauchy stability of (C8a) implies that the roots 2 satisfy 

|ki| < 1, {< 2 ] ^ ^ I ^ i ^ ^ 

The only problem is that 

z ■ 1 when ic = 1 

Therefore, we have to prove whether there is any generalized eigenvalue 
(J. Oliger and B. Gustafsson, private communication) for (C9) and (CIO). For 
stability, we do not want \<i\ 1 from below as | z | -► 1 from above. 

Therefore we want to find out if = 1. Let z*=l-H6, 6>0, and 
K = 1 -1“ £ with 6,e small, we get 

(1 + £)(1 + 6 - 1) = A(1 + 6)[(1 + e)^ - 1] 

^ - Xe(2 + e) 

(1 + e) - Xe(2 + e) 

Since 6 > 0, this implies e > 0; thus, < = ! + £: is There- 

fore z = 1 is not a generalized eigenvalue, and the entire scheme is 
unconditionally stable . 

By applying the same procedure, it can be shown that the boundary approxi- 
mation (C8b) , that is, spatial linear extrapolation, together with the interior 
schemes (a) central difference in space and (b) two-step backward Euler in 
time, form an unconditionally stable scheme for the model equation (Cl). 
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APPENDIX D 


MATRIX METHOD AND POSITIVE REAL FUNCTION METHOD 


Consider a scalar hyperbolic 


equation: 


|^+c|^=0 0<x<l. t>0' 

u(x,0) » f(x) c > 0 
u(0,t) = g(t) 


(Dl) 


The above equation constitutes a well-posed problem. Let vj ■ v(jAx,t) be 
the difference solution of (Dl) at x - jAx, where Ax is the step size. 

Let us discuss the method of lines approach by using central difference in 
space. We will examine the stability of this difference scheme by the matrix 
method (ref. 29) and by the positive real-function method (ref. 30) for fixed 
Ax. 


A word of caution; these methods only show the stability of the ODE for 
a fixed Ax, In order to show that the original PDE is stable, the related 
ODE has to be stable as Ax 0. That is, additional analysis is required. 
The additional requirement involves the analysis of infinite dimension 
matrices. Here, we only show the method for fixed Ax, and want to point out 
that stability of the ODE for fixed Ax does not rule out the possibility 
that the ODE might become unstable as Ax 0. 

By central differencing in space, (Dl) becomes' 



At the right boundary (the numerical boundary) we use the backward difference 
scheme. 


8x 


Ax 


and, therefore, we have 


In matrix notation 



(D2) 


If the real part of all the eigenvalues of A are negative, we can apply a 
stable ODE solver to integrate (D2). The particular type of ODE solver 
depends heavily on the spectrum of the eigenvalues of A, that is, on the 
stiffness of the system. The matrix A cannot be transformed to a diagonally 
dominant matrix with all its diagonal elements positive. We cannot get an 
explicit bound for Ax, We have to actually compute the eigenvalues of A. 

Gary (ref. 29) has shown that A is a stable matrix for various mesh spacings. 


We now turn to the use of positive real functions in an investigation of 
numerical stability of (D2) with fixed Ax. For details of the theory, please 
refer to Dahlquist^s original paper (ref. 30) on this subject. 


Let z * (2Ax/c)X, with X the eigenvalues, N the dimension of A, and 


Dj^(z) = det(zl - A) = 0 


Then D^(z) is of the form 




z i 0 

-1 z 1 

-1 z 1 


-1 z 1 
-2 z + 2 
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It is easy to see that 


D (z) = zD (z) + D (z) n > 3 
n+ 1 n n- i 

If it turns out that D^(z) + 0 for Re(z) > 0 and that imaginary zeros 
are simple, then for each n, all solutions of the ODE’ s are bounded, and any 
A-stable method can be used for the integration in time. 

Let 


D 


n 


D 


n-i 


then 


"n+i 


Z + T- n > 3 


DjCz) 

DjCz) 


z 1 0 

-1 z 1 

0 -2 z + 2 

1 I 


-2 z + 2 


or 


<^3 


z + 2 

^ z(z + 2) + 2 


(D 3 ) 


Since 


(z) »z+2=0 Z--2 

D2 (z) * z(z +2)+2«0 z»-l±i 

have their only zeroes in the left half-plane, it is sufficient (though not 
necessary) to show that are positive functions for n > 3 . Let us look 

at the second part of (D 3 ). Recall that for an arbitrary complex number W 
that if Re(W) > 0 then Re(W”^) > 0 . Since 


f(z) 


z (z + 2) + 2 . 2 

TTl “ ^ ■*'TTT 


is a positive function, it follows that <1)3 is a positive function. By apply- 
ing the proof by induction, we can easily show that n > 3 , are 


64 


positive functions. Thus, the central scheme is stable for c > 0 for fixed 

Ax. 
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APPENDIX E 


SUFFICIENT STABILITY CONDITIONS 


From the discussion of well-posedness of hyperbolic initial boundary /alue 
problems, no new difficulties arise if we have smooth variable coefficients 
and quasi-linear equations with smooth solutions. We will concentrate on a 
general strictly hyperbolic constant coefficient system with well-posed coupled 
boundary conditions. With this system in mind, we will give a detailed 
description in the following order; (1) the basic idea, (2) dissection of the 
problem, (3) reduction of the system to scalar equations, and (4) sufficient 
stability conditions. 


Basic Idea 

The sufficient stability conditions only involve properties of methods 
for related Cauchy problems. We want stable schemes for the related Cauchy 
problem applied at the interior, and stable and uncentered dissipative schemes 
for the related Cauchy problem applied at the boundary. The stabilities of 
the related Cauchy problems are usually known or can be verified by standard 
techniques. The main theories behind these are based on the Cauchy stability 
of the composite method, and the matching of stable schemes, which has been 
examined by Ciment (ref. 13) and Oliger (ref. 15). The usefulness of these 
results is fourfold: (1) stability can be easily verified by standard tech- 

niques; (2) the result can be used to guide us in the construction of stable 
methods for the entire problem; (3) the Cauchy stability of the composite 
method is especially useful and efficient for higher order schemes; and (4) the 
result can help to simplify the verification of the necessary and sufficient 
conditions tremendously if the use of higher order schemes is desired. 


Dissection of the Problem 


We will discuss the approximation of the well-posed strictly hyperbolic 
system 


w + Aw = 0 

t X 

w(x,0) = f(x) 

Liw(O.t) = fi(t) 
L2W(l,t) = f2(t) 


0<x<l, t>ol 


/ 


(El) 


where A is a N N constant matrix, and Li and Lz are rectangular 
matrices. After an appropriate nonsingular transformation T, we can trans- 
form (El) into 
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Ut + = 0 

u(x,0) = f(x) 


(E2a) 


where 



u « Tw 

with boundary conditions 

u^(O.t) = SjU^^(O.t) + gj(t) 

u^^(l,t) = SjjU^(l.t) + gjj(t) 


(E2b) 


where Sj is an «,x(N - 1) matrix, Sn is a (N - £) x («•) matrix, and 
gl = Tf ^ = Tf^, 

From the well-known theorem of Gustafsson et al. (ref. 9), the stability 
of the approximation for an initial boundary value problem on 0 < x < 1 is 
equivalent to the stability of two related quarter-plane problems. Therefore 
we can split (E2) into the related left and right quarter-plane problems. The 
related right quarter-plane problem on 0<x<°°,t>0 is obtained by simply 
removing the boundary at x *■ 1 and extending the definition of our initial 
data and interior approximation to x = that is, 


+ Au^ = 0 

0<x<”, t<0\ 

u(x,0) = f(x) 

1 (E3) 

1^(0, t) = S^u^^(0,t) + gj(t) 

) 
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The related left quarter-plane problem on 


—GO 


<x<l,t>0 is defined 


+ Au^ =0 -®<x<l,t>0^ 

u(x,0) - f(x) 


,E4) 


Therefore, the discussion of the right quarter-plane problem (E3) is suffi- 
cient for our purpose. The analysis of (E4) is similar. 


Reduction of the System to Scalar Equations 

We want to solve system (E3) using finite difference schemes. Divide the 
x-axis into subintervals of length Ax and the t-axis into subintervals of 
length At, Denote the grid points by Xj = jAx and grid functions by 
vj(t) “ v(xj,t), t = nAt and approximate (using one step in time for illus- 
tration; theory holds for multistep) (E3a) in the interior of the domain by 

P 

V (t + At) = 2 j = r, r + 1, r + 2, . . . (£5) 

J i=_p J 

where p is the order of the spatial differencing for the interior scheme, 
and the approximation grid points are defined as in figure L2 (without the 
right boundary present) and are fixed N x n diagonal matrices. 


LEFT 

BOUNDARY 

POINTS 



INTERIOR POINTS 


1 


1 


RIGHT 

BOUNDARY 

POINTS 


0 r-1 r r+1- 


J J+1 J+K 


Figure 12.- Grid point definition 

For the outflow unknowns (variables with negative eigenvalues), we 
approximate the boundary conditions by the following uncentered scheme 





+ At) 


s 

E 

i =-m 




j .m • 0, 


r - 1 , m < j 


(E6) 


where are fixed diagonal (N - £) x (N - 1) matrices and s is the 
order of the spatial differencing for the boundary scheme. Note that for 
m = 0 the scheme is one-sided. 
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The following are a few of the spatially one-sided and uncentered schemes: 

3v^(t) Vj^^(t) - Vj(t) 

3x Ax 

3 Vj(t) - j Vj(t) 

3x Ax 

8Vj(t) + 6vj^^(t) - 3Vj(t) - 2Vj_^(t) 

3x 6Ax 

The first two are one-sided and are of order of accuracy Ax and Ax , respec- 
tively. The last one is uncentered and is of order Ax^. 

For the inflow part (variables with positive eigenvalues), we have the 
analytical boundary condition 

VQ^(t) - SjV^^(t) + ij(t) (E7) 

together with r - 1 additional approximations of the form 

V.^(t) - ^ D vj^(t) + g,(t) (E8) 

J i-o J 


where are fixed £x(N - 1 ) matrices, q is a positive integer, and the 

g^(t) are vectors depending on Ax and |j(t). See Goldberg and Tadmor 
(ref. 14) for derivation of (E8) . 


Since the Aj are diagonal, we can split the scheme (E5) into its inflow 
and outflow parts (ref. 14): 


P 

Vj^(t + At) - j - r, r + 1, . . 

i— p 


Vj^(t + At) 


E 

i w— n 


Ai Vj^^(t) 


j - r, r + 1, . . 


(E9) 


where 
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Now we can see that equations (E5) and (E6) can be partitioned into the 
lowing problems 


(t + At) 


r 

■E 


*1 


j = r. r + 1, . . 

.ElOa) 


v/(t) = SjvJ^t) + |j(t) 


(ElOb) 


Vj^(t) 


T. 

■E 

i=o 


DjiV^^(t) + gj(t) 


J"l» • • -jr—l 


(ElOc) 


vj^(t + ^ j = r, r + 1. . . . 

i=-p 


>(E11) 


E 


i=-m 


° ^ =■ 0, . . .. r - 1, m < j 

i=*“tn 


The outflow problem (Ell) is self contained, while the inflow problem 
(ElO) depends on the outflow part to the extent that the outflow computations 
provide the inhomogeneous boundary values in (ElOb) and (ElOc) . Therefore the 
stability of the right quarter plane under the above approximation is equiva- 
lent to the following two separate parts. 


1. Stability of the inflow problem (ElO) with inhomogeneous boundary 
values 

2. Stability of the outflow problem (Ell) 

Since all the and cj^^ are diagonal matrices, the inflow prob- 

lem splits into I independent approximations and the outflow problem splits 
into (N - «,) independent approximations. Similarly, we can split our left 
quarter-plane problem into the equivalent form. Therefore, the stability 
study of a system of the form (El) reduces to a study of a single scalar equa- 
tion with two related quarter-plane problems as follows 
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+ cti^ - 0 0 S X < ® 

u(x,0) ■ f(x) t > 0» c > 0 

u(0,t) - fej(t) 

+ cu^ - 0 -« < X ^ 1 

u(x,0) * f(x) t i 0, c > 0 

For c > 0, or 

+ cu^ - 0 0 < X < ~ 

u(x,0) ■ f(x) t > 0, c < 0 

+ cu^ *0 -» < x < 1 

u(x,0) - f(x) t > 0. c < 0 

u(l,t) - 

For c < 0. 

Sufficient Stability Conditions 

Let us assiime c > 0» and discuss stability analysis of difference 
approximations to equations (E12). Using the same difference approximation as 
before 

P 

Vj (t + At) ■ ^i^j+i^^^ j • Tf r + 1, . • . 

i— p 

v^(t) - Sj(t) 

Vj(t) * r>l 

P 

Vj(t + At) - j 5 j (E15a) 

i-p 

m ® 

'y ] + At) - . . .,J + K (E15b) 

i—8 j <m + JSJ + K 
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J 


where the approximation grid points are defined as in figure 12 (without the 
left boundary present), and the Ai, Cji are now scalar constants and g.(t) 

^ expansions of the solution in the neighborLod 

of the boundary in terms of the physical boundary data ii(t). This has been 
shown to acquire the desired accuracy of order d if the data is sufficiently 
smooth. The form of gj(t) is as follows; cientiy 


gj(t) 


■E 

i»o 


~^ iT^ —j^u(O.t) + 0(Ax‘*'*'^) 
3x^ 


ij(t) + 0(Ax‘^+^) 
dt 


i“0 

K .-i” Vj(t), j - 1, . . r - 1 can be approximated or extrapolated 

by other uncentered methods (see Oliger, ref, 15). But a stability proof will 
be more complicated. Now, the stability of the inflow, right quarter-plane 
problem (E14) is an immediate consequence of the stability of the Interior 
approximation, so the stability discussion will only deal with the outflow 
left quarter-plane problem (E15) . We need the following definition and 
assumptions: 

Definition; An approximation is said to be Cauchy stable if It 
for the related Cauchy problem. ^ ^ stable 

Assumptions; (1) We assume that our interior approximations and boundary 
approximations are stable for the related Cauchy problems; (2) we assume our ^ 
sl^?ive)?^‘”'°’^”““°''® dissipative (or at least one of the scheme is dis- 

The sufficient conditions rest on the following three results; 

ref 15) theory of matching of stable schemes (Clment, ref. 13; Oliger, 


2. The theory of successive!/ 
posite method (Oliger, ref. 10). 


constructing Cauchy stable methods - com- 


of Gustafsson et al. (ref. 9) - if the method is Cauchy 
stable, then it is stable for the left quarter-plane problem. 

. Matching of stable schemes- If a Cauchy stable scheme of the form (E15a) 
fori (E15M dissipative approximation of the 

st™li ^^Thls li hiJr ^ ^5® resulting approximation is Cauchy 

^ the result of Ciment and Oliger 's theorem on the 

>^he Cauchy stability 

of both methods and the dlssipatlvity of at least one method. The result is ^ 
best illustrated by figure 13. result is 
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I 1 1: 






CAUCHY STABLE 
SCHEME 


CAUCHY STABLE AND 
DISSIPATIVE SCHEME 



COMBINED CAUCHY STABLE SCHEME 


Figure 13.- Sufficient condition. 

Successively constructing Cauchy stable methods - By applying the previous 
method of "matching of stable schemes" on n^ ■ J (see fig. 12), with scheme 
(E15a) for j < J and scheme (El5b) defined for j ■ J -P 1 for all 
j > J 1» the combined scheme is Cauchy stable. We can construct a second 
composite method using the one ve have just constructed with scheme (£l5a) for 
j < J + 1 and the scheme (E15b) defined for j » J + 2 for all j > J + 2. 
This in turn again is Cauchy stable by the method of matching of stable 
schemes. We proceed in this way until we get to j - J + K. This is illus- 
trated by the diagram in figure 14. 

Theory of Gustafsson et al.- By successive construction of a Cauchy stable 
scheme using the composite method, and the assumption we made for (E15) , the 
result of Gustafsson et al. (ref. 9) says that the left quarter-plane (outflow) 
problem is stable. 

Therefore, the key to constructing stable schemes for the initial bound- 
ary value problem for the hyperbolic equations is to have Cauchy stable schemes 
for the interior points and the boundary points, and at least one of the 
schemes is dissipative. 
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J 


J+1 


J+K 


X 


J-t-2 


STEP 1 CAUCHY STABLE ^ 
SCHEME (Eq. E15 (a) ) 


CAUCHY STABLE AND DISSIPATIVE SCHEME 
(Eq. E15(b)) DEFINED FOR] -J+1 


STEP 2 CAUCHY STABLE SCHEME 
OF COMPOSITE METHOD 
FROM STEP 1 


CAUCHY STABLE AND DISSIPATIVE SCHEME 
(Eq.E 15(b)) DEFINED FOR j-J+2 


STEP K+1 


J+K- 1 J+K 


CAUCHY STABLE SCHEME 
OF COMPOSITE METHOD 
FROM STEP K 


CAUCHY STABLE AND 
DISSIPATIVE SCHEME 
(Eq. E15(b)) DEFINED FOR 
i- J+K 


y 

COMBINED CAUCHY STABLE SCHEME 


Figure 14.- Successively constructing Cauchy stable methods. 
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